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Preface 

The  purpose  of  this  study  was  to  investigate  the  optimization  of  stochastic 
response  surfaces  with  linear  programs.  It  was  found  that  using  the  traditional  approach 
of  estimating  a  response  surface  and  using  it  as  the  objective  function  of  a  linear  program 
yielded  a  bias  in  the  mean  solution.  Also,  nonoptimal  extreme  points  have  a  large 
probability  of  b'^ing  chosen.  This  research  investigated  a  method  to  overcome  these 
disadvantages  and  obtain  an  improved  solution. 
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Abstract 


V 

This  rese-irch  investigated  an  alternative  to  the  traditional  approaches  of 
optimizing  a  stochastic  response  surface  subject  to  constraints.  This  research 
investigated  the  bias  in  the  expected  value  of  the  solution,  possible  alternative  decision 
variable  settings,  and  a  method  to  improve  the  solution.  A  three  step  process  is  presented 
to  evaluate  stochastic  response  surfaces  subject  to  constraints.  Step  1  is  to  use  a 
traditional  approach  to  estimate  the  response  surface  and  a  covariance  matrix  through 
regression.  Step  2  samples  the  objective  function  of  the  linear  program  (i.e.,  response 
surface)  and  identifies  the  extreme  points  visited.  Step  3  presents  a  method  to  estimate 
the  optimal  extreme  point  and  present  that  information  to  a  decision  maker. 


1  Introduction 


1.1  Background 

Constrained  and  unconstrained  optimization  is  a  diverse  and  growing  field  with 
applications  in  many  areas.  Two  key  optimization  areas  of  interest  in  this  lesearch  were; 
optimization  of  simulation  output  with  constraints  and  stochastic  programming. 

Response  surface  methodology  can  be  used  to  optimize  a  system  modeled  by  a 
simulation.  Response  surface  methodology  typically  incorporates  three  areas: 
experimental  design,  regression  analysis,  and  optimization.  The  response  surface  (i.c., 
the  fitted  regression  model)  generated  from  a  stochastic  simulation  is  a  metamodel  that 
approximates  the  system  being  simulated.  The  coefficients  of  the  response  surface 
representing  the  simulation  in  a  design  region  are  random  variables.  The  random,  or 
stochastic,  nature  of  the  coefficients  in  the  response  surface  is  at  the  core  of  this  research. 
Typically,  after  establishing  and  validating  a  response  surface  the  analyst  ignores  its 
stochastic  nature  and  employs  it  to  estimate  optimum  operational  conditions  in  a  resource 
constrained  environment  (11:138). 

Stochastic  programming,  usually  a  completely  separate  field  from  response 
surface  methodology,  concentrates  on  the  stochastic  nature  of  elements  in  math 
programming  problems.  One  aspect  of  stochastic  programming  concentrates  on  the 
random  nature  of  the  coefficients  in  the  objective  function  in  a  linear  program  while 
assuming  the  constraints  are  deterministic.  A  literature  review  suggests  there  is  little 
incorporation  of  stochastic  programming  to  simulation  optimization;  the  research  done  is 
limited  in  scope  and  solves  only  individual  problems.  No  research  effort  has  been  found 
that  investigates  the  general  process.  Davis  and  West  observe  that: 
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...recent  research  has  demonstrated  that  there  is  indeed  a  role  for  the  employment 
of  mathematical  programming  procedures  with  simulations...  in  modeling 
situations  where  unceilainties  are  small,  mathematical  programming  would  be  the 
prefen'ed  choice.  However,  many  decisionmaking  problems  do  not  meet  this 
criterion.  For  example,  long-term  production/inventory  planning  is  often 
analyzed  using  a  mathematical  programming  formulation  even  though  this 
problem  is  known  for  its  uncertainties  arising  in  long-term  demand  forecasts,  in 
the  costs  of  input  material,  and  in  manufacturing  productivities...  Another  classic 
problem  for  which  considerable  uncertainties  often  exist  is  the  investment 
portfolio  selection  problem...  For  any  problem  with  severe  uncertainties  the 
adopted  approach  of  merging  simulation  with  mathematical  programming  can  be 
applied.  (5:200,209). 


1.2  Problem  Statement 


The  goals  of  this  research  arc  threefold.  First,  this  research  will  characterize  the 
impact  of  the  stochastic  nature  of  the  response  surface  on  the  optimization  of  a 
simulation  with  multiple  constraints.  Second,  this  research  will  develop  a  means  to 
identify  the  "true"  optimum  point  in  practice.  Third,  this  research  will  develop  a  method 
to  increase  confidence  in  the  point  estimate  on  the  optimal  answer  to  the  constrained 
optimization  problem. 

1.3  Objectives  &  Related  Methods 

The  goals  of  this  research  was  met  by  accomplishing  the  following  objectives. 
Objective  1:  To  characterize  the  impact  of  the  stochastic  nature  of  the  objective  function 
(response  surface)  coefficients  on  the  distribution  of  (the  optimal  solution  to  the  linear 
programming  problem)  and  identify  both  the  basis  and  extreme  point  changes. 


Methods: 

1.  Generated  a  variety  of  linear  programming  problems  to  investigate  the  topic. 

2.  Wrote  a  computer  program  that  performs  a  Monte  Carlo  process  on  the  linear 
program  by  varying  the  objective  function  coefficients  and  error  term  to  be 
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sampled.  Collected  data  on  the  distribution  of  basis  changes,  arid  sampled 
optimum  extreme  points. 

3.  Defined  a  variety  of  variance-covariance  matrices  (a2(X'X)'l)  to  generate  the 
multivariate  normal  distribution  of  the  coefficients. 

4.  Defined  five  noise  levels  to  test. 

5.  Ran  a  Monte  Carlo  procedure  and  evaluate  the  mean  and  variance  for  each 
extreme  point,  basis,  and  for  the  overall  solution. 

Objective  2:  Investigate  a  method  to  design  an  experiment  on  the  response  surface 
coefficients  that  samples  the  response  stirface  in  an  effective  way. 

Method:  Used  a  Box-Bchnken  design  and  investigated  modifications  to  the  design  that 
better  sampled  the  response  surface. 

Objective  3:  Evaluated  ways  to  better  estimate  the  "true"  optimal  extreme  point,  and 
present  the  data,  to  a  decision  maker,  in  an  efficient  way. 

Mctliods: 

1 .  Used  a  ranking  and  selection  procedure  to  screen  extreme  points  and  identify 
the  "best"  extreme  point. 

2.  Investigated  presentation  of  data  through  a  histogram. 


3 


2  Literature  Review 


2.1  Introduction 

The  following  paragraphs  will  review  literature  pertinent  to  optimization  of  a 
stochastic  response  surface  subject  to  multiple  constraints.  Specifically,  the  discussion 
covers  the  topics  of  experimental  design,  constrained  simulation  optimization,  and 
stochastic  programming  with  Monte  Carlo  simulation. 

2.2  Discussion 

2.2.1  Experimental  Design 

Much  research  has  been  done  on  experimental  designs  since  the  1950s  (7:571). 
The  basic  goal  of  experimental  design  is  to  choose  the  settings  of  the  factor  levels  in  a  set 
of  experiments.  In  simulation,  experimental  design  is  the  process  of  developing  a 
scheme  to  conduct  an  experiment  on  a  simulation  and  collect  output  information  deriving 
the  maximum  amount  of  useful  information  with  a  minimum  expenditure  of  resources. 

In  experimental  design  jargon,  the  input  variables  to  the  simulation  are  factors  and  the 
levels  (values)  those  variables  can  take  are  treatments. 

Experimental  design  in  simulation  is  distinguished  from  experimental  design  in 
general  by  two  things.  First,  in  simulation  the  analyst  specifies  the  factors  and  treatments 
at  the  beginning  to  get  the  optimal  design.  Second,  the  input  random  variables  for  each 
simulation  experiment  are  controlled  by  the  analyst  who  can  exploit  this  by  making 
comparisons  between  experiments  more  precise,  such  as  using  the  same  pseudo-random 
number  stream  for  different  experiments  (i.e.,  using  common  random  numbers). 

Typical  objectives  of  experimental  design  in  a  simulation  are  to 

1.  Understand  the  effect  of  the  factors  on  the  experimental  output. 

2.  Estimate  the  parameters  of  interest. 
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3.  Make  a  selection  from  a  set  of  alternatives. 

4.  Find  the  treatment  levels  for  all  factors  that  produce  the  optimum  response 
(3:105-223). 

There  arc  many  types  of  designs,  from  the  simple  2^  factorial  to  specialty  designs 
to  measure  lack  of  fit,  fit  second  order  models,  or  control  variance.  The  list  of 
experimental  designs  is  extensive,  although  "the  central  composite  [design]  is  used  more 
than  any  other  family  of  RSM  designs"  (11:147). 

A  key  influence  in  the  stability  in  predicted  variance  is  the  rotatability  of  the 
design:  while  a  high  degree  of  rotatability  is  desirable,  a  perfectly  rotatability  design  is 
not  needed  (11:138).  Rotatability  refers  to  the  way  variance  propagates  through  the 
design;  in  a  perfectly  rotatable  design  variance  is  only  a  function  of  the  distance  from  the 
center  of  the  design.  In  other  words,  in  a  perfectly  rotatable  design,  variance  can  be 
defined  as  increasing  in  concentric  circles  expanding  out  from  the  center  of  the  design. 

In  a  near-rotatable  design  the  variance  can  be  thought  of  in  the  same  way  except  the 
circles  arc  not  perfect,  but  slightly  elliptical.  In  light  of  variance  considerations,  Myers, 
Khuri  and  Carter  conclude 

...the  RSM  user  needs  to  learn  from  the  Taguchi  approach  that  system  variability 
should  be  a  major  component  in  the  analysis.  A  similar  argument  can  be  made 
for  consideration  of  the  distribution  of  variance  of  the  prediction  in  the 
assessment  of  experimental  designs.  Often  the  success  of  the  RSM  endeavor  is 
dependent  on  the  properties  of  y  [hat]  at  different  locations  in  the  design  space. 
Many  standard  designs  have  prediction  variances  which  increase  dramatically  as 
one  gets  close  to  the  design  parameter.  As  a  result,  any  conclusions  drawn 
(regarding  choice  of  optimal  or  improvements  in  operating  conditions) 
concerning  response  near  the  design  boundary  are  suspect.  Yet  we  see  very  little 
that  deals  with  this  in  design  assessment  or  comparisons  among  designs.  We  too 
often  evaluate  a  design  on  the  basis  of  one  number  (say,  D-efficiency)  when  the 
important  aspects  of  behavior  are  multidimensional.  (1 1 : 152) 


Biles  and  Swain  observe  "the  n-dimensional  simplex  design,  which  employs  n+1 
design  points  at  the  vertices  of  a  regular  simplex,  gives  the  greatest  efficiency  in  terms  of 
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information  per  design  point  "  (2:138).  Experimental  design  is  a  diverse  area;  only 
through  investigation  will  identification  of  the  family  of  "best"  designs  occur. 


2.2.2  Constrained  Simulation  Optimization 

Classical  simulation  optimization  is  a  broad  field.  Typical  approaches  to 

optimization  include  heuristic  searches,  complete  enumeration,  random  searches,  steepest 

ascent,  coordinate  searches,  pattern  searches  and  many  othcis  (6:117-121). 

Biles  and  Swain  present  several  strategics  for  constrained  simulation  optimization 

that  appear  to  be  representative  of  the  main  constrained  simulation  optimization  efforts 

going  on  today.  They  fit  and  validate  a  response  surface  using  an  n-dimcnsional  simplex, 

biradial,  or  cquiradial  design.  They  account  for  the  variance  of  the  error  term,  but  they 

assume  the  "response  surfaces  are  the  expected  values  of  the  observed  responses." 

(2:135).  They  do  not  account  for  the  stochastic  nature  of  the  response  surface,  but 

employ  a  recursive  method  by  applying  an  optimization  procedure  and  then  returning  to 

the  simulation  model  until  the  optimal  criteria  are  met.  Their  procedures  include  direct 

search  techniques,  first-order  response  surface,  and  second-order  response  surface 

procedures.  The  type  of  constraints  Biles  and  Swain  used  are  either  simple  (upper  and 

lower  bounds)  or  multiple  (e.g.,  budget  or  resource)  constraints  (2:135-137). 

The  choice  between  point  estimation  verses  interval  estimation,  while  important 

to  the  analyst,  has  only  recently  been  addressed: 

Many  users  of  RSM  allow  conclusions  to  be  drawn  concerning  the  nature  of  a 
response  surface  and  the  location  of  optimal  response  without  taking  into  account 
the  distributional  properties  of  the  estimated  attributes  of  the  underlying  response 
surface.  Although  the  distribution  of  these  quantities  has  not  been  considered 
directly,  efforts  have  been  made  to  develop  interval  estimates.  Box  and  Hunter 
(1954)  used  a  version  of  Fieller's  theorem  to  develop  a  100(1  -  alpha)% 
confidence  region  for  the  location  of  the  stationary  point.  The  construction  of  a 
confidence  interval  around  the  response  at  the  stationary  point  of  the  true  surface 
has  only  recently  been  a  subject  of  interest  in  the  statistical  literature.  Khuri  and 
Conlon  (1981)  gave  an  expression  for  the  bounds  of  an  interval  conditional  on  the 
estimated  location  of  the  stationary  point.  (11:146) 
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The  experimental  design  impacts  the  distributional  form  of  the  underlying  data,  which  in 
turn,  "has  an  impact  on  the  estimation  of  the  model  parameters  and  on  the  inferences 
drawn  from  an  RSM  analysis"  (11:146). 

Constrained  optimization  of  stochastic  simulations  through  linear  programming  is 
starting  to  appear  in  the  literature,  but  no  literature  has  been  found  that  evaluates  the 
impact  of  the  stochastic  nature  of  the  response  surface  in  the  general  case.  Myers,  Khuri 
and  Carter  support  this:  "Nearly  all  practical  RSM  problems  are  tnily  multiple  response 
in  nature.  Sophisticated  ways  of  solving  [stochastic]  multiple-response  problems  are  not 
generally  well  known,  however"  (11:147). 

2.2.3  Stochastic  Programming 

Stochastic  programming  considers  three  stochastic  areas  of  math  programming: 
the  objective  function  c^x,  constraint  matrix  A,  and  the  right  hand  side  b  vector. 
Traditionally,  stochastic  programming  has  not  been  incorporated  in  simulation 
optimization  analysis. 

There  are  three  classical  approaches  to  solving  stochastic  programming  problems 
where  the  coefficients  vary:  expected  value,  "Fat,"  and  "Slack"  approaches.  The 
expected  value  app.roach  uses  the  point  estimate  for  the  coefficients  to  solve  the  math 
programming  problem.  The  "Fat"  approach  chooses  a  pessimistic  value  for  the 
coefficients  in  the  math  programming  problem.  The  "Slack"  solution  method  assumes 
randomness  in  the  constraint  matrix  and  the  right  hand  side  and  adds  a  penalty  function 
to  the  objective  function.  The  "Slack"  method  is  a  two-stage  problem  that  assumes  the 
decisionmaker  can  adjust  a  previous  decision  (9:463-470).  These  three  approaches  to 
sta'hastic  programming  are  used  when  the  coefficients  are  either  random  or  constrained 
to  a  given  set.  Another  area  of  stochastic  programming  is  chance  constrained 
optimization:  "[In  this]  approach,  it  is  not  required  that  the  constraints  should  always 
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hold,  but  we  shall  be  satisfied  if  they  hold  in  a  given  proportion  of  cases  or,  to  put  it 
differently,  if  they  hold  with  given  probabilities"  (13:75). 

Most  analytical  methods  assume  no  basis  changes  take  place,  and  many  do  not 
account  for  changes  in  the  extreme  point  solutions  (12:21 1),  Under  these  assumptions 
analytical  methods  have  been  developed  to  estimate  the  mean  and  the  variance  of  the 
stcichastic  program: 

...[Consider]  the  cost  of  minimizing  C  =  c'^'x  [the  objective  function],  where  the 
components  of  the  vector  c  had  a  joint  normal  distribution  with  means  mj  and 

covariance  matrix  V.  C  is  then  also  normally  distributed  with  mean  M  =  m^x 
and  variance  =  x^Vx. 

If  all  coefficients  and  constants  of  the  consti'aints  are  fixed,  and  if  we  want  to 
minimize  the  expected  value  of  C,  then  the  problem  is  reduced  to  the 

deterministic  program  of  minimizing  m^^x.  But  let  us  assume  that  we  want  to 
minimize  the  expected  utility  of  C,  which  we  define  as:  1  -  exp(  -aC)  where  a  is  a 
positive  constant,  which  economists  call  a  measure  of  the  aversion  to  risk. 
(13:23-24) 

Vajda  goes  through  a  proof  that  reduces  the  problem  to  a  "deterministic  program  of 
minimizing  aM  -  .5a2s2,  a  function  which  is  quadratic  in  x"  (13:24). 

The  feasibility  of  the  solution  to  a  stochastic  program  is  a  topic  of  much  research. 
Randomness  in  the  constraint  matrix  or  the  right  hand  side  b  vector  can  cause  the 
solution  to  be  super-optimal  and  infeasible.  If  only  the  objective  function  is  stochastic, 
the  solution  will  remain  feasible  and  only  the  optimal  value  will  change  (13:3). 

Bard  and  Chatterjee  introduce  a  perspective  of  the  variability  of  a  design  that  is 
specifically  concerning  objective  function  bounds  for  the  inexact  linear  programming 
problem  with  generalized  cost  coefficients.  They  conclude  that  increasing  the 
"oblongness"  of  the  variance  of  the  coefficients  in  the  objective  function  rather  than 
decreasing  their  volun.e  provides  better  results  when  solving  the  stochastic  linear 
program  (1:491).  In  other  words,  uniform  variance  in  the  coefficients  of  the  objective 
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function  is  not  necessarily  desirable,  but  decreasing  variance  of  one  coefficient  at  the 

expense  of  others  may  lead  to  improved  LP  solutions. 

Bracken  and  Soland  present  a  paper  on  "a  statistical  decision  analysis  of  a  one- 

stage  linear  programming  problem  with  deterministic  constraints  and  stochastic  criterion 

function"  (4:205).  This  paper  presents  analytical  and  Monte  Carlo  methods  to  find  the 

expected  value  of  both  perfect  and  sampled  information,  but  it  is  not  possible  to  solve  the 

analytical  problem  when  the  objective  function  coefficients  come  from  a  multivariate 

normal  distribution.  In  addition,  this  article  presents  a  method  to  describe  the 

"...distribution  of  the  optimal  value  of  the  linear  programming  problem  with  stochastic 

objection  function  and  [discusses]  Monte  Carlo  and  numerical  integration  procedures  for 

estimating  [the  optimal  value]  distribution"  (4:205). 

A  misconception  exists  that  there  is  no  need  to  evaluate  the  stochastic  nature  of 

the  objective  function  in  a  linear  program  because  the  analyst  can  use  sensitivity  analysis 

to  conduct  a  proper  evaluation.  Davis  and  West  address  this  misconception: 

Post-optimal  or  sensitivity  analysis  provides  the  modeler  with  the  ability  to 
analyze  the  functional  behavior  of  the  optimal  solution  as  parameter  assignments 
are  modified,  but  these  methods  again  provide  little  insight  toward  the  probable 
values  that  the  optimal  solution  [of  a  stochastic  problem]  will  assume.  For  this 
reason  it  is  often  difficult  to  choose  a  robust  solution.  (5:199) 

Bracken  and  Soland  obsei’ve  that  the  optimal  solution  to  any  linear  program  will 
always  occur  at  an  extreme  point:  therefore  only  extreme  points  need  to  be  considered. 

In  other  words,  if  all  the  extreme  points  can  be  identified,  then  the  actual  linear 
programming  problem  docs  not  need  to  be  solved,  but  instead  a  Monte  Carlo  procedure 
can  be  used  to  evaluate  the  probability  distribution  of  the  simulation  output  at  each 
extreme  point.  The  value  of  the  extreme  points  can  then  be  evaluated  and  compared  to 
identify  the  top-ranking  alternatives.  This  article  refers  to  an  article  by  C.  E.  Clark  in 
which  he  describes  a  method  (using  a  reduced  set  of  extreme  points)  approximating  the 
characteristics  of  the  distribution  of  the  "...maximum  value  of  the  (reduced)  linear 
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programming  problem  admitting  only  the  p  selected  alternatives."  (4:212,214).  "p  is  a 
subset  of  the  overall  (assumed  independent)  extreme  points. 

There  are  several  reasons  for  studying  the  reduced  linear  programming  problem. 
First  of  all,  the  number  r  of  extreme  points  which  comprise  the  set  S  may  be  so 
large  that  the  decision  maker  finds  it  undesirable  to  consider  ail  of  them  in  his 
analysis  of  the  decision  situation.  He  may  therefore  decide  to  limit  his  further 
analysis  to  the  p  selected  alternatives.  The  smaller  number  of  alternatives  thus 
available  makes  computation  of  the  EVPI  [Expected  value  of  perfect  information] 
and/or  the  EVSI  [Expected  value  of  sampled  information]  much  more  feasible. 
Second,  the  expected  value  of  the  maximum  value  of  the  reduced  problem  is  a 
lower  bound  to  the  conesponding  quantity  for  the  original  problem,  and  could 
therefore  serve  as  an  estimate  of  the  expected  value  (>f  the  maximum  value  of  the 
original  problem.  We  would  expect  the  bound  to  be  best  when  the  p  selected 
alternatives  are  the  p  top-ranking  alternatives  (in  order)  with  respect  to  the  prior 
distribution  on  c  [the  unknown  mean  vector  of  a  multivariate  stochastic  process]. 
(4:212) 


Clark's  procedure,  as  contrasted  with  a  Monte  Carlo  simulation 

has  the  advantage  that  the  results  may  be  obtained  quickly  and  cheaply  for 
different  selections  of  the  alternative  extreme  points  and/or  different  parameters 
in  tlie  distribution  of  c.  The  accuracy  achieved  with  Clark's  procedure  is 
somewhat  limited,  however,  especially  when  the  distribution  of  v  [where  v  is  the 
objective  function]  is  degenerate,  whereas  great  accuracy  can  be  achieved  in 
Monte  Carlo  simulation  if  a  sufficiently  large  number  of  draws  is  uscd...Onc 
difficulty  is  that  the  values  obtained  vvith  Clark's  procedure  are  dependent  upon 
the  order  in  which  the  variables  are  listed.  (4:214-224) 


Davis  and  West  present  both  a  decision  theory  approach  and  a  Monte  Carlo 
approach  to  solving  this  problem. 

Decision  theory  begins  by  predefining  the  alternative  solutions  that  will  be 
considered.  Next  potential  realizations  for  the  decisionmaking  environment  must 
be  specified  with  apriori  assignment  of  the  probability  that  each  realization  occur. 
Using  this  information,  the  tiade-offs  among  the  alternatives  are  then  analyzed, 
and  the  apparent  optimal  solution  alternative  is  selected.  Decision  theory  does 
generate  the  probability  that  each  alternative  will  be  the  optimal  solution. 
However,  the  analysis  is  limited  by  the  number  of  predefined  alternative  solutions 
selected  and  the  accuracy  in  the  specification  of  the  a  probabilities  for  the 
potential  states  of  the  systems...  generation  of  probabilistic  bounds  upon  the 
optimum  solution  requires  considerable  computational  effort  with  the 
mathematical  programming  approach  whereas  decision-theory  approaches 
provide  this  information  directly...Although  these  analyses  are  still  typically 
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limited  to  the  consideration  of  a  finite  set  of  alternatives,  the  methods  of  response 
surface  methodology  coupled  with  classical  optimization  approaches  have  been 
used  to  optimally  assign  key  parameter  values  within  simulated  systems  operating 
under  the  selected  alternative.  (5:199-200) 


Davis  and  West,  instead  of  defining  the  possible  solution  states,  use  the  simulated 
decisionmaking  scenarios  from  the  Monte  Carlo  analysis  and  assume  each  scenario  has 
an  equal  probability  of  occurring  (5:207).  The  above  approach  belies  the  difficulty  in 
identifying  the  possible  scenarios  and  their  probability  of  occurring.  The  paper  was 
introduced  as  a  contrast  between  the  methods  of  decision-theory  and  Monte  Carlo 
simulation,  but  the  decision  analysis  builds  on  the  results  from  the  Monte  Carlo  analysis. 
Also,  the  paper  presents  an  analysis  of  only  a  single  study  and  not  an  evaluation  of  the 
overall  process  (5:199-209). 

2.3  Summary 


Incorporating  stochastic  programming  in  the  optimization  procedure  for 
stochastic  simulations  v,ith  constraints  has  recently  been  used  in  limited  cases,  but  no 
overall  evaluation  of  the  process  has  been  donek  In  solving  problems  with  random 
coefficients 

difficulties  arise  from  two  sources.  First,  meaningful  simulation  models  must  be 
generated  from  which  the  uncertain  parameter  values  can  be  sampled.  Second, 
the  ability  to  statistically  analyze  the  results  from  numerous  sample  optimal 
solutions  generated  during  the  simulation  must  be  demonstrated.  Even  if  both 
difficulties  are  addressed,  the  task  still  remains  of  selecting  the  decision  which 
provides  the  best  compromise  between  optimality  and  risk.  (5:209) 


Morben,  in  solving  a  "real  world"  problem,  demonstrates  a  case  where  using  the  expected 
value  of  a  stochastic  objective  function  leads  to  an  answer  which  falls  outside  a  95% 
coi  'idence  bound  found  through  a  Monte  Carlo  analysis  (10:27).  This  case  clearly 


'  T5is  research  will  analyizc  the  overall  process  and  look  into  ways  to  incorporate  experimental  design  in 
the  stcchastic  programming  subportion. 
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demonstrates  there  is  a  risk  in  some  situations  if  only  the  expected  value  is  used,  and  it 
makes  the  case  for  incorporating  some  form  of  stochastic  analysis.  Experimental  design, 
constrained  simulation  optimization,  and  stochastic  programming  arc  all  well  developed 
fields  while  the  need  to  incorporate  aspects  of  all  three  seems  to  exist  little  research  has 
been  done  to  achieve  this  goal. 
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3  Phase  I:  Impact  of  Estimation  Errors  Methodology 

3.1  Introduction 


Phase  I  investigated  how  errors  in  estimating  the  response  surface  affect  the 
estimate  of  the  optimal  solution  of  the  maximizing  linear  program.  Figure  1.  shows  the 
basic  approach 


This  research  phase  considers  the  simulation  a  black  box  that  consists  of  a  'Tnith 
Model"  plus  noise. 
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Figure  2.  Black  Box  Simulation  Model 


The  simulation  response  serves  as  input  to  the  estimate  of  a  response  surface  (the 
objective  function  of  a  linear  program).  The  Functions: 

Z*  =  LP(C,A,b)  (1) 

:^*  =  LP(t:,A,b)  (2) 

define  the  optimal  value  TT  (or  estimated  optimal  value  of  a  linear  program 
employing  the  revised  simplex  method.  Where 

C  =  c^x  true  (or  known)  objective  function 
t  =  c^x  estimated  objective  function 
A  =  constraint  matrix 
b  =  right  hand  side  vector. 

c  =  true  surface  coefficients  underlying  the  metamodel 
c  =  c  +  8  estimated  coefficients  of  objective  function  (response  surface) 
Assumption:  E~N(0,a2(xTx)'l) 

Phase  I  started  with  the  premise  that  7^  (the  "true"  optimum)  is  equal  to  the 
expected  value  of  the  parameters  of  the  linear  program: 

Z*  =  LP(E(^,A,b))  (3) 

but,  it  will  be  demonstrated  that  7^  is  not,  in  general,  equal  to  the  expected  value  of  the 
linear  program  with  the  estimated  value  of  the  objective  function: 

Z*  =  LP(C,  A,  b)  =  LP(E(^,  A,  b))  ^  E(LP(^,  A,  b))  =  E(^*)  (4) 
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Funher,  as  the  standard  error  in  the  estimates  of  the  coefficients  increases  the  bias  in 
and  increases. 

3.3  Testing  Hypothesis 

This  research  used  a  computer  program  to  test  the  hypothesis  in  Equation  4.  The 
computer  program  uses  a  Monte  Carlo  approach;  sampling  from  a  "Truth  Model"  with 
noise  it  generates  a  response  surface  which  is  used  as  the  objective  function  of  a  linear 
program.  The  objective  function  is  then  sampled  using  the  variance-covariance  matrix 
generated  from  the  regression  of  the  design  matrix  and  response  while  collecting 
statistics  at  all  stages  (see  Figure  3). 

3.4  Investigating  Indicated  v.  True 

Besides  characterizing  the  distribution  of  tlie  estimated  optimal  solution,  this 
research  investigated  what  solution  an  analyst  might  expect  in  practice  versus  the  "true" 
solution.  The  analyst  estimates  two  key  elements:  the  optimal  extreme  point  (T^*)  and 
the  optimal  value  given  that  point  (^*).  Theoretically,  any  changes  to  the  linear  program 
could  change  characteristics  of  the  comparison,  but  this  research  will  investigate  whether 
there  is  a  common  trend  to  be  identified. 

To  characterize  the  possible  results  the  computer  program  generates  an  estimated 
objective  function  (response  surface)  and  solves  the  linear  program.  Both  the  estimated 
optimal  value  and  the  estimated  optimal  extreme  point  are  generated  and  compared  to  the 
true  optimal  extreme  point,  its  value,  and  the  value  of  the  tine  function  at  the  estimated 
extreme  point  (later  these  values  are  sorted  and  plotted  to  give  a  visual  representation  of 
the  comparison).  Contrasting  these  plots  with  plots  of  different  linear  programs  and  at 
different  noise  levels  gives  insights  to  the  problem. 
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4  Results  Phase  I 


4.1  Introduction 

There  is  no  singular  "true"  answer  to  Phase  I,  but  it  appears  the  general  process 
can  be  characterized.  Consistent  results  were  found  for  all  maximization  problems 
studied  by  evaluating  different  linear  programs  at  different  noise  levels  and  generating 
random  "Truth  Models".  The  following  results  capture  what  can  be  expected  in  the 
general  case.  As  noise  is  introduced  in  the  estimation  of  the  objective  function  the 
estimated  optimal  extreme  point  will  vaiy  as  shown  in  Figure  4. 


An  additional  module  was  added  to  the  computer  program  to  generate  random 
"truth  models,"  constraint  matrices,  and  right  hand  side  vectors.  The  output  from  this 
module  characterized  this  basic  approach  over  a  variety  of  linear  programming  problems. 
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Always,  for  the  randomly  generated  maximization  problems  tested,  the  computer 
program  indicated  that  the  value  of  the  "True"  linear  program  is  less  than  or  equal  to  the 
value  of  the  linear  program  with  noise : 

Z*  <  E(LP(^,  A,  b))  =  E(i*)  (5) 

This  can  be  graphically  represented  in  the  following  figure: 


Figure  5.  Noise  Impact  &  Bias 


this  figure  shows  how  a  normally  distributed  estimate  of  the  objective  function  t  affects 
the  normal  distribution  in  the  estimate  of  the  optimal  value 

The  high  bias  in  the  estimate  of  the  mean  was  present  in  all  linear  programming 
problems  analyzed  (and  those  randomly  generated).  Also,  as  the  noise  level  increased  in 
a  given  problem  the  bias  increased  in  a  roughly  linear  trend;  here  the  bias  is  the  mean 
estimated  optimum  minus  the  true  optimum. 

Bias  =  E(2*)  -  Z*  (6) 

The  following  figures  illustrate  a  typical  case  where  the  standard  en'or  =  a(parameter 
estimates).  The  underlying  linear  program,  in  this  case,  is: 
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Bias 


Figure  6.  Bias  Inflation 


The  bias  increases  as  the  standard  error  increases,  and  the  standard  deviation  of 
follows  a  similar  trend  as  shown  in  the  Figure  7. 


Combining  the  last  two  points,  one  can  expect  as  a(parameter  estimates)  increases  not 
only  does  the  bias  increase,  but  the  spread  in  :^'^increases. 

4.3  indicated  v.  Actual  Solutions 

The  bias  in  the  estimated  value  of  the  optimal  answer  only  illustrates  the  trend 
over  multiple  realizations,  but  an  analyst  may  be  interested  in  what  he  can  expect  in  one 
realization  of  the  process.  Here,  let's  go  a  step  beyond  the  statistics  in  the  previous 
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section  dealing  with  bias  and  analyze  the  contrast  between  each  and  the  actual  value 
of  choosing  that  extreme  point.  In  the  following  figures 

True  =  value  of  "Truth  Model"  given  the  extreme  point  con’esponding  to 
Indicated  = 

The  extreme  point  for  is  that  extreme  point  the  linear  program  chooses  as  optimal. 

The  following  figures  are  characteristic  of  all  problems  evaluated.  The  following  points 
are  of  interest: 

1.  The  True  solution  is  sorted  in  descending  order  where  the  first  (left  most) 
extreme  point  represents  the  "true"  optimal  extreme  point. 

2.  The  Indicated  solutions  are  sorted  in  descending  order  around  their  extreme 
point. 

3.  All  points  in  the  plots  are  equally  likely  in  a  single  realization  of  the  process 
(8000  samples). 

4.  For  a  point  estimate  there  is  no  way  to  know  if  the  estimate  is  high  or  low 
without  prior  knowledge. 

5.  The  distribution  of  the  solution  around  each  extreme  point  is  normally 
distributed  with  a  bias. 

6.  As  o  increases  the  percent  of  visits  to  the  "hue"  optimal  extreme  point 
decreases. 

7.  As  a  increases  there  is  a  greater  chance  that  the  chosen  extreme  point  will  be 
greatly  inferior  to  the  "true"  optimal  extreme  point,  but  there  is  still  a  chance 
with  be  much  higher  than  even  the  "true"  optimal  extreme  point  Z^. 

8.  The  chance  of  selecting  the  solution  with  the  chosen  extreme  point  within  10% 
using  occurs  only  about  10%  of  the  time. 

9.  Extreme  point  changes  occur  where  a  True  value  step  change  occurs  and 
where  the  Indicated  value  jumps  from  a  low  value  to  a  high  value. 
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When  a  =  2.25  visits  to  the  "true"  optimal  extreme  point  occur  about  14%  of  the 
time,  and  about  98%  of  the  solutions  are  "close"  to  the  "true"  optimal  extreme  point. 


Figure  9  Indicated  Z*  vs  Actual  a  =  3.25 
When  a  =  3.25  selecting  the  "true"  optimal  extreme  point  occurs  about  7%  of  the 
time,  and  about  92%  of  the  solutions  are  "close"  to  the  "tnie"  optimal  extreme  point. 
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Figure  10  Indicated  Z*  vs  Actual  a  =  4.25 
When  a  =  4.25  selecting  the  "true"  optimal  extreme  point  occurs  about  6%  of  the 
time,  and  about  91%  of  the  solutions  are  "close"  to  the  "tme"  optimal  extreme  point. 

A  single  realization  of  the  process  (or  point  estimate)  gives  little  information  on 
how  the  chosen  extreme  point  will  actually  perform.  Other  techniques  must  be 
investigated  to  put  this  information  into  practical  application. 
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5  Phase  II  Sampling  Extreme  Points  Methodolo 


5.1  Introduction 

Phase  I  illustrated  that  limited  information  about  the  true  solution  can  be  obtained 
from  a  single  realization  of  the  process.  The  estimated  extreme  point  may  lead  to  a 
highly  biased  solution  compared  to  the  true  extreme  point.  Phase  II  investigated  how  to 
obtain  the  true  extreme  point  using  two  methods.  The  first  method  samples  the  generated 
objective  function  (in  a  Monte  Carlo  fashion)  using  the  variance-covariance  matrix  from 
the  regression  analysis  and  catalogs  the  extreme  points  visited.  The  second  method 
samples  the  generated  objective  function  via  a  design  and  catalogs  the  extreme  points 
visited.  Identifying  the  optimal  extreme  point  may  be  possible  by  sampling  the 
simulation  at  each  extreme  point  visited-this  will  be  investigated  in  Phase  III.  Also, 
empolying  a  screening  technique  will  be  investigated  to  improve  the  efficiency  of  the 
process. 

5.2  Starting  Hypothesis 

As  a  starting  hypothesis,  this  fcscarch  assumes  that  given  an  initial  response 
surface  and  its  associated  variance-covariance  matrix  we  can  sample  the  "true"  optimal 
extreme  point. 

Pr("true"  optimal  EP  sampled  I  initial  estimated  EP  &  G^(X'rX)'l)  »  1 

as  N  -4  oo 

where 

EP  =  extreme  point 

initial  estimated  EP  =  extreme  point  identified  when  the  objective  function  is 
assumed  to  be  deterministic  and  the  LP  is  solved  once. 

N  =  number  of  objective  function  samples. 


o2(XTx)-l  =  variance-covariance  matrix  driving  multivariate  normal  sampling. 

5.3  Visits  to  true  Extreme  Points  via  Sampling 

In  practice,  analysts  usually  only  have  one  estimate  of  the  objective  function,  but 
in  this  research  the  computer  program  can  generate  any  number  of  objective  functions. 
Calculating  a  variance-covariance  matrix  for  each  response  surface  the  program  then  uses 
it  as  input  to  generate  a  multivariate  normal  sample.  The  computer  program  then 
samples  from  this  distribution  an  arbitrary  number  of  times  and  tests  if  the  computer 
program  samples  the  "true"  extreme  point.  The  computer  program  evaluates  this  process 
by  combining  these  two  sampling  routines  in  a  Monte  Carlo  fashion  (call  this  the  brute 
force  approach).  Of  course,  without  knowing  the  "true"  extreme  point  a  priori  there  is  no 
way  of  knowing  the  "true"  extreme  point  was  sampled  in  practice.  In  this  research  the 
"true"  solution  is  known  at  the  outset  because  the  computer  program  defines  the 
underlying  truth  model. 

5.4  Visits  to  true  Extreme  Points  via  Design 

As  an  alternative  to  the  bmte  force  approach  of  Monte  Carlo  sampling  this 
research  investigated  sampling  using  a  design.  This  is  investigated  by  applying  a  Box- 
Behnken  design  (3:519)  with  only  one  sample  at  the  zero  level.  The  goal  is  to  sample  the 
true  extreme  point,  a  design  is  used  to  try  to  minimize  the  number  of  samples  of  the 
objective  function  of  the  linear  program.  Modifications  to  the  Box-Behnken  design  were 
investigated. 

5.5  Screening  Extreme  Points 

Evaluating  the  estimated  objective  function  using  either  Monte  Carlo  sampling  or 
a  design,  requires  solving  a  linear  program  for  every  sample  of  the  objective  function. 
Solving  multiple  linear  programs  requires  a  lot  of  computer  time,  and  for  a  large  problem 


threatens  to  make  this  research  impractical.  If  a  way  could  be  found  to  test  each  sampled 
objective  function  to  see  if  its  optimum  basis  has  previously  been  sampled  then  a 
screening  method  could  be  applied  to  reduce  the  number  of  linear  program  solutions 
required.  In  this  research,  only  the  objective  function  is  stochastic  and  therefore  only 
optimality,  and  not  feasibility,  is  an  issue.  The  optimality  condition,  for  a  maximization 
problem,  in  the  general  case  is: 

c-CbB-1a<0  (7) 

where 

c  =  objective  function 

Cg  =  coefficient  of  the  basic  variables 

B'l  =  inverse  of  columns  under  basic  variables 
A  =  constraint  matrix 
The  specific  equation  is: 

t-tijB-lA<0  (8) 

where 

tf  =  the  coefficients  of  the  sampled  objective  function 

Cg  =  the  coefficients  of  the  sampled  objective  function  associated  with  the  basic 

variables. 

As  new  extreme  points  are  sampled  their  corresponding  B'l  A  matrix  and  a  vector 
identifying  the  basic  variables  are  stored  in  a  set,  this  set  is  then  used  to  screen  new 
objective  function  samples.  For  every  new  sample  of  the  objective  function  the  computer 
program  cycles  through  the  B'l  A  matrices  using  Equation  8  until  the  optimality 
condition  is  satisfied,  if  the  optimality  condition  is  not  satisfied  it  solves  the  linear 
program  to  identify  a  new  basis.  Using  this  technique  a  linear  program  is  solved  only 
once  for  each  unique  extreme  point  sampled. 
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6.1  Introduction 


Sampling  the  estimated  objective  function  to  include  the  "true"  optimal  extreme 
point  in  the  sampled  set  proved  an  effective  method.  While  it  is  not  possible  to  identify 
the  "true"  optimal  extreme  point  when  sampled,  it  is  possible  to  sample  a  set  of  extreme 
points  which  include  the  "tnie"  optimal  extreme  point. 

6.2  "True”  Extreme  PointVisits  w/  Monte  Carlo  Sampling 

Monte  Carlo  sampling  proved  a  effective  method  of  sampling  over  a  broad  range 
of  problems.  As  is  expected,  the  probability  of  sampling  the  "tine"  optimal  extreme 
point  increases  as  the  number  of  Monte  Carlo  samples  increases.  The  numbers  of  Monte 
Carlo  samples  required  for  a  given  confidence  level  presumably  increases  as  the  number 
of  decisions  variables  increases  and  as  the  en'or  in  the  estimate  increases.  This  research 
has  not  investigated  any  method  to  determine  a  prbri  the  number  of  samples  required  to 
achieve  a  given  confidence  level.  The  following  tables  are  typical  results  for  problems 
involving  4  decision  variables.  The  percent  "true"  optimal  extreme  point  miss  is 
evaluated  by  replicating  the  process  1000  times. 


Standard  En'or 

.25 

1.25 

2.25 

3.25 

4.25 

100  Samples 

%  Miss  "tnie" 

.8 

1.8 

2.6 

5.2 

10.0 

200  Samples 

%  Miss  "true" 

1.0 

1.6 

1.9 

2.4 

4.25 

300  Samples 

%  Miss  "true" 

.8 

1.3 

1.5 

2.3 

3.7 

500  Samples 

%  Miss  "true" 

.5 

.9 

1.3 

1.9 

3.1 

Table  1.  Monte  Carlo  Samples 
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An  advantage  to  Monte  Carlo  sampling  is  that  if  the  analyst  has  the  time  and  resources, 
and  wants  to  be  conservative,  the  number  of  Monte  Carlo  samples  could  be  increased. 
While  Monte  Carlo  sampling  is  the  least  efficient,  if  the  analyst  is  willing  to  take  enough 
samples,  it  could  be  the  most  effective  in  sampling  the  "true"  optimal  extreme  point. 

6.3  "True”  Extreme  Point  Visits  w/  Design  Sampling 

First,  a  modified  Box-Behnken  design  was  used  to  sample  the  objective  function 
of  the  linear  program.  The  modification  involved  only  one  sample  at  the  zero  level.  The 
sampling  was  done  by  varying  the  estimated  objective  function  coefficients  by  a  percent 
of  the  estimated  standard  deviation  (called  standard  deviation  multiplier)  in  a  method 
prescribed  by  the  design.  This  approach  was  possible,  in  this  case,  because  an  onliogonal 
design  was  used  to  sample  the  original  "Black  Box  Simulation"  to  estimate  the  response 
surface,  as  a  result  there  are  no  off-diagonal  elements  in  the  variance-covariance  matrix. 
A  more  complicated  method  is  need  if  off-diagonal  elements  were  present,  but  it  seems 
an  initial  orthogonal  design  is  a  reasonable  approach.  An  example  of  a  four-variable  Box- 
Bchnkcn  design  is  found  in  Appendix  B. 

Table  2  illustrates  the  results  from  a  single  Box-Behnken  design  with  shown 
standard  deviation  multipliers. 


Standard  Error 

.25 

1.25 

2.25 

3.25 

4.25 

Standard  Dev 

1.5 

%  Miss  "true" 

2.0 

3.2 

6.7 

18.2 

30.7 

Standard  Dev 

mm 

%  Miss  "time" 

.9 

1.3 

IS) 

22.3 

37.1 

Standard  Dev 

2.5 

%  Miss  "time" 

.2 

.8 

11.5 

32.7 

41.6 

Standard  Dev 

3.0 

%  Miss  "true" 

.1 

1.5 

18.4 

38.9 

50.6 

Table  2.  Single  Box-Behnken  Design 
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Tests  using  this  single  Box-Behnken  design  showed  limited  success.  It  appears  this 
single  design  is  inadequate  to  sample  the  "true"  extreme  point. 


The  second  modification  to  the  standard  Box-Behnken  design  was  to  double  the 
length  of  the  design  by  sampling  at  each  design  point  twice.  For  every  identical  pair  of 
design  points  different  standard  deviation  multipliers  were  used.  In  effect,  a  three-level 
design  was  transformed  into  a  pseudo  five-level  design.  It  is  not  a  true  five-level  design 
because  each  design  point  has  only  three  levels,  it  is  really  the  same  design  run  twice 
with  two  different  standard  deviation  multipliers.  Table  3  shows  results  of  this 
technique. 


Standard  Error 

.25 

1.25 

2.25 

3.25 

4.25 

Standard  Dev's 

m 

%  Miss  "true" 

4.0 

5.5 

6.5 

7.5 

20.5 

Standard  Dev's 

1.5,2 

%  Miss  "true" 

.9 

1.5 

3.6 

9.9 

20.2 

Standard  Dev's 

1.5, 2.5 

%  Miss  "true" 

.2 

.1 

2.8 

10.3 

19.6 

Standard  Dev's 

1.5,3 

%  Miss  "true" 

.25 

0.0 

1.65 

9.1 

20.3 

Standard  Dev's 

1,2.5 

%  Miss  "true" 

.2 

.1 

2.4 

9.9 

17.1 

Standard  Dev's 

1,2 

%  Miss  "true" 

.9 

1.5 

3.6 

8.1 

16.0 

Standard  Dev's 

1,3 

%  Miss  "true" 

0.0 

0.0 

2.2 

11.0 

20.3 

Table  3.  Double  Box-Behnken  Type  Design 


The  double  Box-Behnken  design  (sampling  the  objective  function  49  times) 
shows  promise.  Results  with  the  double  Box-Behnken  design  are  superior  to  sampling  in 
a  Monte  Carlo  fashion  49  times. 


Standard  Error 

.25 

1.25 

2.25 

3.25 

1  4.25 

%  Miss  "tme" 

1.2 

4.0 

5.8 

12.3 

18.6 

Table  4.  49  Monte  Carlo  Samples 


Results  over  a  broad  range  of  problems  indicate  this  design  is  superior,  but  not 
dramatically,  to  an  equivalent  number  of  Monte  Carlo  samples.  In  general,  either  case 
fails  to  give  confidence  in  the  results. 

The  next  modification  includes  adding  a  third  Box-Bchnken  design  to  the 
previous  two  designs  and  sampling  it  at  a  different  standard  deviation,  this  is  a  pseudo 
seven-level  design.  In  essence,  this  is  equivalent  to  sampling  from  three  consecutive 
designs. 


Standard  Enor 

.25 

1.25 

2.25 

3.25 

4.25 

Standard  Dev's 

.5, 1.5, 

2.5 

%  Miss  "true" 

0.3 

0.2 

1.4 

5.8 

10.6 

Standard  Dev's 

1.0, 2.0, 

3.0 

%  Miss  "true" 

0.1 

0.2 

0.8 

4.8 

11.5 

Standard  Dev's 

1.0, 1.75 

2.5 

%  Miss  "tme" 

0.3 

0.2 

1.1 

4.4 

9.9 

Standard  Dev's 

.5, 1.75, 

3.0 

%  Miss  "tme" 

0.1 

0.2 

0.8 

5.0 

10.6 

Standard  Dev's 

1.0, 1.75 

2.5 

°/o  Miss  "tme" 

0.9 

3.2 

7.4 

15.5 

26.1 

Table  5.  Triple  Box-Behnken  Type  Design 


The  triple  Box-Behnken  design  had  good  results,  but  required  more  samples.  In 
this  case,  the  triple  Box-Behnken  design  (with  four  decision  variables)  was  sampled  73 
times.  As  a  comparison,  the  results  of  73  Monte  Carlo  samples  are  presented  in  Table  6. 


Standard  Error 

.25 

1.25 

2.25 

3.25 

4.25 

%  Miss  "tme" 

1.0 

2.2 

4.3 

warn 

Table  6.  73  Monte  Carlo  Samples 
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Again  there  is  an  advantage  to  the  design  over  the  equivalent  number  of  Monte 
Carlo  samples.  The  main  advantage  to  a  Monte  Carlo  approach  is  that  the  number  of 
samples  can  be  arbitrarily  increased  to  achieve  the  confidence  desired,  this  may  be 
desirable  if  a  higher  confidence  in  the  solution  in  needed  than  is  possible  with  this 
design.  To  this  point,  each  design  was  an  improvement  over  an  equivalent  number  of 
Monte  Carlo  samples,  but  no  design  gave  a  high  success  rate  at  higher  noise  levels. 

In  an  effort  to  improve  the  success  rates  with  higher  levels  of  noise  another  type 
of  modification  to  the  basic  Box-Bchnken  design  was  investigated.  In  this  case,  the  basic 
structure  at  each  design  point  was  modified.  Instead  of  sampling  at  the  design  points 
using  a  thrcc-lcvcl  approach  of  1,  -1 ,  or  0,  this  new  design  was  a  true  five-level  design 
where  each  design  point  was  sampled  with  some  combination  of  1 ,  -1 ,  .5,  -.5,  or  0.  This 
modification  doubles  the  length  of  the  design  and  at  each  design  point  alternatively 
samples  form  cither  1  or  .5.  An  example  of  this  new  design  is  found  in  Appendix  B. 


Standard  Error 

.25 

1.25 

2.25 

3.25 

4.25 

Standard  Dev 

1,5 

%  Miss  "true" 

5.1 

6.8 

10.3 

12.5 

15.5 

Standard  Dev 

2.0 

%  Miss  "true" 

1.7 

3.1 

4.0 

6.3 

10.8 

Standard  Dev 

2.5 

%  Miss  "trtic" 

.5 

1.8 

2.1 

5.3 

11.2 

Standard  Dev 

3.0 

%  Miss  "trtic" 

0.0 

.6 

1.7 

E9Hi 

12.3 

Table  7.  Single  5-level  Box-Behnken  Type  Design 


The  single  modified  5-lcvcl  Box-Bchnken  design  has  49  design  point,  the  same 
number  as  the  double  Box-Behnken  design  presented  in  Table  3.  The  5-levcl  design  has 
a  higher  success  rate  in  sampling  the  "tnie"  optimal  extreme  point  than  either  the  double 
Box-Behnken  design,  or  a  equivalent  number  of  Monte  Carlo  samples.  The  5 -level  Box- 
Bchnken  design  represents  an  improvement  when  sampling  at  higher  noise  levels,  but  the 
errors  could  still  be  considered  significant. 
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A  further  modification  attempts  to  decrease  the  errors  in  sampling  the  "true" 
optimal  extreme  point  by  doubling  the  design  and  choosing  a  different  standard  deviation 
multiplier  for  the  second  half  of  the  design.  This  modification  is  analogous  to  the  change 
creating  the  double  Box-Behnken  design.  This  design  creates  a  pseudo  nine-level  design. 
The  results  of  1000  replications  of  this  design  are  contained  in  Table  8. 


Standard  Error  .25  1.25 

2.25  3.25  4.25 

Standard  Dev  .5, 1.75 

%  Miss  "tnic" 

2.67 

4.0 

6.6 

9.2 

12.9 

Standard  Dev  1.0, 2.0 

%  Miss  "tnic" 

1.8 

2.9 

4.5 

6.2 

10.2 

Standard  Dev  1.5, 2.5 

%  Miss  "tme" 

0.8 

1.2 

1.9 

3.2 

6.1 

Standard  Dev  1.5, 3.0 

%  Miss  "true" 

0.2 

0.3  1.0 

2.6  6.8 

Table  8.  Double  5-level  Box-Behnken  Type  Design 


The  double  modified  5-level  Box-Behnken  design  gave  excellent  results.  This  design 
gave  the  best  results  for  methods  with  about  97  samples,  and  it  is  competitive  with  a 
Monte  Carlo  method  of  200  samples. 

In  the  next  modification  another  modified  5-level  design  is  added  and  sampled  at 
a  different  standard  deviation.  This  pseudo  13-lcvel  design  (four  variables)  has  145 
design  points.  The  results  are  found  in  Table  9  and  show  excellent  results. 


Standard  EiTor 

.25 

1.25 

2.25 

3.25 

4.25 

Standard  Dev 

1.0, 2.0, 

3.0 

%  Miss  "true" 

0.1 

0.3 

0.3 

1.9 

3.0 

Standard  Dev 

1.5, 2.5, 

3.5 

%  Miss  "tme" 

0.0 

0.1 

0.1 

1.1 

2.9 

Standard  Dev 

1.5, 2.75 

4.0 

%  Miss  "tme" 

0.0 

0.0 

0.1 

0.5 

1.9 

Standard  Dev 

1.5, 3.0, 

4.5 

%  Miss  "tme" 

0.0 

0.0 

0.0 

0.5 

2.2 

Table  9.  Triple  5-level  Box-Behnken  Type  Design 
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145  samples  of  the  triple  5 -level  Box-Behnken  design  was  superior  to  all  Other  designs 
and  even  superior  to  500  Monte  Carlo  samples.  This  design  provides  excellent  sampling 
in  an  relatively  efficient  manner.  The  main  drawback  is  that  it  requires  145  samples  with 
only  four  variables.  If  this  full  design  is  run  its  length  will  double  with  every  added 
variable.  As  a  possible  way  to  offset  the  time  required  to  solve  this  many  linear 
programs  this  research  investigated  a  method  of  screening  sampled  objective  functions  to 
decrease  the  computations  required. 

Another  interesting  consideration  is  the  number  of  extreme  points  visited  with 
different  sampling  techniques.  If  one  sampling  method  provided  high  accuracy,  but 
required  may  more  extreme  points  to  be  sampled,  then  it  might  not  be  the  best  design  to 
employ.  Fortunately,  no  design  greatly  increased  the  number  of  extreme  points  sampled. 
Table  10  illustrates  the  total  unique  extreme  points  sampled  for  200  Monte  Carlo  samples 
and  two  design,  these  results  arc  typical  of  all  sampling  options. 


Standard  Eiior 

.25 

1.25 

2.25 

3.25 

4.25 

200  Monte  Carlo 

#  unique  ext.  points 

3 

4 

6 

8 

9 

Double  5-level  Box- 
Behnken  Type 

1.5, 2.5 

#  unique  ext.  points 

2 

3 

5 

8 

8 

Triple  5-level  Box- 
Behnken  Type 

1.5, 2.75, 

4.0 

2 

4 

5 

7 

8 

Table  10.  Total  Unique  Extreme  Points  Sampled  by  Case 


6.4  Screening  Extreme  Points 

The  main  drawback  to  using  the  above  approaches  is  the  need  to  solve  a  linear 
program  for  every  sampled  objective  function.  Screening  the  new  objective  function 
samples  proved  an  efficient  technique.  With  this  technique,  a  linear  program  is  solved 
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only  once  for  each  sampled  extreme  point.  The  improved  efficiency  will  vary  from 
problem  to  problem  and  will  also  depend  on  the  number  of  objective  function  samples, 
but  improvement  can  be  measured  in  orders  of  magnitude.  Applying  this  technique 
greatly  increases  the  practicality  and  efficiency  of  the  research.  Using  this  screening 
procedure  makes  a  strong  case  for  using  the  triple  5-level  Box-Behnken  design  approach 
with  a  large  number  of  samples. 

In  this  research  the  IMSL  Fortran  Library  was  used  to  evaluate  the  revised 
simplex  method,  but  it  does  not  return  a  B’l  matrix,  only  values  for  the  optimal  value 
(Z"^),  the  primal  and  dual  solutions  of  the  decision  variables.  To  employ  this  screening 
technique  the  matrix  must  be  found.  The  B  matrix  was  found  by  choosing  the 
columns  under  the  basic  variables.  The  values  of  slack  variables  are  not  given;  therefore, 
when  slack  variables  are  basic  the  computer  program  uses  the  principle  of 
complcmentaiy  slackness  to  identify  the  values  of  the  coiiesponding  columns  in  the  B 
matrix. 
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7  Phase  III:  Selecting  the  Optimal  Extreme  Point 

7.1  introduction 

Phase  II  established  a  method  to  sample  extreme  points  with  the  goal  of  including 
the  "true"  extreme  point  in  the  sample.  Phase  Ill's  aim  is  to  pick  the  "tme"  optimal 
extreme  point  from  the  population  of  sampled  extreme  points.  To  do  this,  consider  the 
extreme  points  as  settings  of  the  input  decision  variables  to  the  simulation. 

7.2  Starting  Hypothesis 

After  identifying  the  feasible  extreme  points,  the  linear  program  is  no  longer 
needed,  and  the  decision  variable  settings  at  any  extreme  point  are  used  as  input  to  the 
simulation  to  estimate  Z*.  Once  a  set  of  decision  variable  settings  is  selected  the 
problem  becomes  selecting  the  "best"  option.  This  research  starts  with  the  premise  that 
given  a  set  of  sampled  extreme  points  containing  the  "true"  optimal  extreme  point,  it  can 
be  identified  as  the  optimal  given  enough  sampling: 

<I)i  =  E(EPi)  as  N  -4  oo  (9) 

IfEP*cEPthcn 

E(EP*)  =  max  (hj  =  Z*  (10) 

i 

where 

N  =  number  of  samples  from  simulation  (i.e.,  simulation  inns) 

Oj  =  Value  of  'Truth  Model"  at  ith  extreme  point 

EP'''  =  "true"  optimal  extreme  point 

EP  =  set  of  sampled  extreme  points 

EPj  =  ith  extreme  point  in  sample  (arbitrary  ordering) 

Using  this  method  it  may  be  possible  to  both  identify  the  "tme"  optimal  settings  for  the 
decision  variables  (extreme  point)  and  an  unbiased  estimate  for  the  optimal  solution  7^. 
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But,  this  method  requires  many  samples  from  the  simulation  and  only  considers  the 
means  and  not  the  distributions  around  the  means. 

7.3  Ranking  and  Selection  of  Decision  Variable  Sets 

Law  and  Kelton  present  ranking  and  selection  procedures  that  offers  an 
alternative  to  the  brute  force  method  presented  above  (8:596).  The  analyst  may  be 
interested  in  three  selection  approaches.  First,  the  procedure  is  based  on  selecting  the 
best  of  the  k  decision  variable  settings.  Second,  as  an  initial  screening  procedure, 
selecting  a  subset  of  size  m  that  contains  the  best  of  the  k  decision  variable  settings. 
Third,  selecting  the  best  m  alternatives  from  the  k  decision  variable  settings,  this 
approach  would  offer  greater  flexibility  for  the  decision  maker  by  providing  more 
options. 

7.3.1  Selecting  the  Best  of  k  Systems 

Previously,  represented  the  estimated  optimal  solution  from  the  (inear pw^mnu 
Now,  let  the  Zy's  represent  point  estimates  from  the  simufaiion  for  the  jth  replication  of 

the  ith  decision  variable  set  and  pj  =  E(Zjj).  This  approach  assumes  the  Zy's  are 

independent  (8:596).  Law  and  Kelton  define  a  method  of  finding  the  smallest  expected 
response;  this  research  focuses  on  the  largest  expected  response  defining  pj^to  be  the  /th 

largest  of  the  pj's  and 

Pij>Pi2>-*->Pij,  (11) 

we  want  P(CS)>P*  (12) 

provided  pij  -  Pi2  ^  d"^  (13) 

and  the  minimal  CS  probability  is  P"^  >  1/k  (14) 

where 

CS  =  correct  solution 

d"^  >  0  defines  the  "indifference"  specified  by  the  analyst. 

Law  and  Kelton  present  a  procedure  originally  developed  by  Dudewicz  and  Dalai  (1975) 
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[Involving]  "two-stage"  sampling  from  each  of  the  k  systems.  In  the  first  stage 
we  make  a  fixed  number  of  replications  of  each  system,  then  use  the  resulting 
variance  estimates  to  determine  how  many  more  replications  from  each  system 
are  necessary  in  a  second  stage  of  sampling  in  order  to  reach  a  decision.  It  must 
be  assumed  that  the  [Zjj's]  arc  normally  distributed,  but  (importantly)  we  need  not 

2 

assume  that  the  values  of  o.  =  Var(Zij)  are  known;  nor  do  we  have  to  assume  that 


2 

the  a.  are  the  same  for  different  i's.  (Assuming  known  or  equal  variances  is  very 

unrealistic  when  simulating  real  systems.)  The  procedure's  performance  should  be 
robust  to  departures  from  the  normality  assumption,  especially  if  the  [Zjj's]  are 

averages.  (8:596) 


Law  and  Kclton  define  the  first-stage  sampling  with  hq  >  2  replications  of  each  k 
decision  variable  sets  and  define: 


for  i  =  1,  2, , . , ,  k.  N  is  total  sample  size  needed  for  system  i 


(15) 


(16) 


f  rh?s?(«o)ii 

Nj  =  max'!  «o+ 1,1  — ;7T~  I  (l'^) 

I  I  IJ 

where  [ x1  is  the  smallest  integer  greater  than  or  equal  to  the  real  number  x,  and  h  j 
(which  depends  on  k,  and  «o)  is  a  constant  that  can  be  obtained  from  Table  1 
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(Appendix  C).  Next,  make  Nj  -  hq  more  replications  of  system  i  (i  =  1,2, . . .  ,k)  and 
solve  for  the  second-stage  sample  means  by 


(Ni-V 


Ni 

E  Zij 

j=«0  +  1 


Ni-no 


(18) 


Define  the  weights  as 


(19) 


and  Wi2  =  1  -  Wji,  for  i  =  1,2, . . .  ,k.  Also,  define  the  weighted  sample  means  as 


Zi(Ni)  -  Wii  z/‘\«o)  +  Wi2  Zi^CNj  -  ng)  (20) 

and  select  the  decision  variable  setting  (extreme  point)  with  the  largest  Zi(Ni)  (8:597). 
Law  and  Kelton  conclude: 

The  choice  of  and  d*  depend  on  the  analyst's  goals  and  the  particular  system 
under  study;  specifying  them  might  be  tempered  by  the  computing  cost  of 
obtaining  a  large  Nj  associated  with  a  large  or  small  d"^.  However,  choosing 
710  is  more  troublesome,  and  we  can  only  say,  on  the  basis  of  our  experiments  and 
various  statements  in  the  literature,  that  itQ  be  at  least  20,  (8:  597-598) 

7.3.2  Selecting  a  Subset  m  Containing  the  Best  of  k  Decision  Variables 

There  may  be  cases  where  EP  (the  set  of  sampled  extreme  points)!  is  large  and 
the  above  approach  would  require  too  much  computer  time  and  effort;  here  an  initial 
screening  procedure  could  be  useful. 
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Now  consider  selecting  a  subset  of  size  m  containing  the  best  from  k  decision 
variable  settings  (i.e.,  EPj  where  i  =  1, 2, . . . ,  k)  where  m  is  specified  by  the  analyst. 

This  is  that  same  situation  as  above  except  "correct  selection  (CS)  is  defined  to  mean  that 
the  subset  of  size  m  that  is  selected  contains  a  system  [extreme  point]  with  mean  {ijj  and 

we  want  P(CS)  >  P"^  provided  that  -  pj2  ^  d"^;  here  we  must  have  1  <  m  <  k  -1,  P"^  > 

m/k,  and  d*>0  (8:599)."  Also,  replace  hi  by  h2  that  depends  not  only  on  k,  P*,  and  «o, 

but  also  on  m  (see  Table  2  Appendix  C). 

Then  we  make  Nj  -  no  niore  replications,  from  the  second  stage  means 

-(2) 

Zj  (Nj  -  Hq),  weights  Wil  and  Wi2,  and  weighted  sample  means  Zi(Nj), 

[exactly  as  before].  Finally,  we  define  tiie  selected  subset  to  consist  of  the  m 

systems  corresponding  to  the  m  smallest  values  of  the  Zj(Ni)'s.  (8:599) 

As  m  increases  "considerably  fewer  replications"  are  required  than  when  m  =  1  (8:600). 

7.3.2  Selecting  the  Best  m  of  k  Decision  Variable  Settings 

Using  the  m  best  of  k  decision  variable  settings  approach,  the  analyst  can  provide 

the  best  m  alternatives  to  the  decision  maker  giving  him  a  broader  base  for  a  decision. 

This  approach  is  very  similar  to  the  above  two  approaches  except  the  subset  of  size  m 
equals  the  largest  expected  responses  pjp  Pi2, . . (this  set  is  not  ranked  ordered). 

Here  P(CS)  >  P*  provided  ^  Also,  P"^  >  m!(k-m)!/k!  and  replace  h2  by 

h3  (see  Table  3  Appendix  C). 

7.4  Histogram  Comparison 

The  ranking  and  selection  procedure  above  presents  a  method  of  evaluating  the 
expected  value  of  different  decision  variable  settings,  but  choosing  the  "best"  solution 
often  involves  more  than  just  identifying  the  largest  expected  value. 

Histograms  are  plotted  using  all  simulation  samples  from  the  best  m  alternatives. 
The  histogram  can  the  aid  the  decision  maker  by  visually  representing  the  possible 
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realizations  of  the  process  at  given  settings.  Two  important  advantages  are:  avoiding  risk 
by  choosing  the  smallest  variance,  and  illustrating  nearly  equivalent  alternatives  and 
allowing  the  decision  maker  to  consider  factors  not  captured  by  the  model.  For  instance, 
figure  1 1  illustrates  possible  histograms  from  the  top  three  decision  variable  sets. 


Figure  11  Histogram  Comparison 


It  is  not  clear  which  is  the  best  alternative.  The  top  plot  has  the  highest  mean,  but  a  risk 
averse  decision  maker  may  choose  the  second  option  to  avoid  the  possible  down  side  of 
the  first  option.  In  either  case,  a  visual  representation  presents  the  decision  maker  with  a 
broader  knowledge  base  from  which  to  make  a  decision. 
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8  Phase  III  Investigation  Results 


A  ranking  and  selection  procedure  was  used  to  analyze  the  extreme  points 
sampled  by  the  double  Box-Behnken  type  design  when  the  standard  error  equals  3.25. 
The  sampled  extreme  points  are  in  Table  11,  and  include  the  actual  at  each  extreme 
point  which,  of  course,  would  not  be  known  in  practice. 


Extreme  Point 

xl 

x2 

x3 

x4 

True  Z* 

1 

4 

0 

0.666667 

6.6666667 

215.333 

2 

1 

11 

0 

0 

212 

3 

2 

10 

0 

0 

210 

4 

2 

8 

0 

2 

216 

5 

8 

0 

2 

0 

166 

6 

0 

8 

2 

0 

182 

7 

0 

0 

4.666667 

2.666667 

147.333 

8 

0 

0 

6 

0 

118 

Table  11 .  Unique  Extreme  Points  Sampled 


Because  there  are  only  eight  extreme  points  no  screening  is  needed,  and  the 

ranking  and  selection  of  the  best  m  out  of  k  alternatives  was  used.  Lets  assume  the  best 
three  alternatives  are  desired  at  a  90%  confidence  level  (i.e.,  h3  =  3.532)  with  d=  4. 

Where  d  is  the  minimum  separation  between  extreme  point  sample  means  desired,  and 
like  the  confidence  level  is  chosen  by  the  analyst.  Knov/ing  the  tine  separation  between 
Z*  it  is  clear  d  is  too  high,  here  it  is  chosen  a  little  high  to  illustrate  the  robustness  of  the 

process.  There  will  always  be  a  trade  off  between  the  confidence  level  desired,  the  value 
of  d,  and  the  number  of  samples  (Nj)  required  at  each  extreme  point.  As  the  confidence 

level  increases,  or  the  value  of  d  decreases,  the  number  of  samples  required  increases. 
Table  12  correctly  identifies  the  top  three  alternatives,  recall  that  the  order  of  the  three 
alternatives  is  not  guaranteed.  In  this  case,  the  decision  maker  would  choose  between 
extreme  points  one,  two  and  four. 
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1 

Zil(20) 

■  Si(20) 

Ni 

Zi2(Ni-20) 

Wii 

Wi2 

Zi(Ni) 

1 

217.876 

154.707 

121 

217.135 

.186 

.814 

217.273 

2 

206.517 

118.576 

93 

212.157 

.247 

.753 

210.767 

3 

209.244 

122.273 

96 

208.107 

.242 

.758 

208.837 

4 

214.886 

254.121 

199 

215.379 

.12 

.88 

215.32 

5 

167.00 

182.121 

142 

164.912 

.142 

.858 

165.209 

6 

186.924 

112.299 

88 

179.201 

.257 

.743 

181.186 

7 

145.86 

165.889 

130 

146.701 

.18 

.82 

146.55 

8 

113.123 

135.033 

106 

116.445 

.221 

.779 

115.791 

Table  12.  Selecting  the  Three  (w/  d=4)  Best  of  the  Eight  Extreme  Points 
Tabic  13  shows  a  case  where  a  different  approach  is  taken.  The  initial  estimates  of  the 

means  Zj^(20)  show  there  appears  to  be  a  natural  division  between  the  first  four  extreme 

points  and  the  last  four.  Because  of  this  natural  division  it  might  be  advantageous  to 

choose  the  best  four  extreme  points  out  of  the  set  of  eight.  In  this  case,  at  a  90% 
confidence  level  h3  equals  3.571  and  d  is  again  chosen  equal  to  four. 


i 

Zil(20) 

S-(20) 

Ni 

Zi2(Ni-20) 

Wii 

Wi2 

Zi(Ni) 

1 

217.876 

154.707 

124 

216.534 

.189 

.811 

216.788 

206.517 

118.576 

95 

212.193 

.24 

.76 

209.244 

122.273 

98 

208.554 

.234 

.766 

208.716 

214.886 

254.121 

203 

215.129 

.133 

.887 

215.102 

167.00 

182.121 

146 

164.80 

.163 

.837 

165.159 

186.924 

112.299 

90 

179.182 

.253 

.747 

181.143 

7 

145.86 

165.889 

134 

146.854 

.178 

.822 

146.677 

113.123 

135.033 

108 

116.624 

.208 

.792 

115.895 

Table  13.  Selecting  the  Four  (w/  cl=4)  Best  of  the  Eight  Extreme  Points 
Table  13  coiTectly  identifies  to  top  four  alternatives.  Even  overestimating  d,  this  method 


proved  effective.  Table  14  illustrates  another  approach,  the  Zil(20)  seem  to  show  two 

distinct  groups.  This  can  be  exploited,  instead  of  d=4  let  us  choose  d=6,  this  is  done  with 

the  goal  of  separating  the  groups  and  then  getting  a  feel  for  the  rankings.  This  approach 
might  be  taken  if  multiple  replications  are  difficult  to  make  (note  the  decrease  in  Nj  from 
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Table  13  to  Table  14).  Also,  this  approach  illustrates  another  advantage  of  this 
technique,  using  the  ranking  and  selection  procedure  it  is  possible  to  identify  the  top  ni 
competing  alternatives.  In  other  words,  this  helps  identify  the  number  of  roughly 
equivalent  alternatives  (i.e.,  the  value  for  m). 


i 

Zi'(20) 

Si(20) 

Ni 

Zi2(Ni-20) 

Wii 

Wi2 

Zi(Ni) 

1 

217.876 

154.707 

215.661 

.393 

.607 

216.531 

2 

206.517 

118.576 

213.92 

.542 

.458 

209.908 

3 

209.244 

122.273 

44 

208.454 

.517 

.483 

208.863 

4 

214.886 

254.121 

91 

215.833 

.263 

.737 

215.583 

5 

167.00 

182.121 

65 

162.854 

.348 

.652 

164.296 

6 

186.924 

112.299 

40 

179.464 

.537 

.463 

183.472 

7 

145.86 

165.889 

59 

148.435 

.369 

.631 

147.484 

8 

113.123 

135.033 

48 

115.913 

.446 

.554 

114.669 

Table  14.  Selecting  the  Four  (w/  d=6)  Best  of  the  Eight  Extreme  Points 


After  choosing  a  ranking  and  selection  method,  the  confidence  level,  and  d  then  this 
research  recommends  using  the  data  obtained  through  the  procedure  to  create  a 
histogram.  The  histogram  is  a  way  to  aid  the  decision  maker.  In  this  example,  all  the 
actual  variances  arc  equal,  but  this  method  has  its  strengths  when  the  variances  arc 
different.  Figure  12  illustrates  the  histograms  of  the  top  four  alteniatives,  the  dotted 
vertical  lines  represent  Zi(Ni)  for  each  alternative.  An  alternative  to  presenting  a 

histogram  of  the  data  is  the  plot  the  nomal  probability  curve  defined  by  the  estimated 
mean  and  varianec. 
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Figure  12  Sample  Case  Histogram  Comparison 


At  this  point  the  choice  is  up  to  the  decision  maker. 
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9  Conclusions  and  Recommendations 


9.1  Introduction 

This  research  presents  a  superior  method  to  the  traditional  approach  of  estimating 
a  response  surface  and  then  using  it  as  the  objective  function  of  a  linear  program.  Over 
multiple  realizations  the  traditional  approach  will  overestimate  the  true  mean  response, 
and  it  is  unlikely  that  the  "true"  optimal  extreme  point  will  be  chosen.  Small  variance  in 
the  estimates  of  the  response  surface  coefficients  can  lead  to  large  variance  in  the 
estimation  of  and  a  low  probability  of  choosing  the  correct  optimal  extreme  point 
EP*.  By  using  the  screening  procedure  this  general  procedure  may  become  practical  for 
general  application. 

9.2  Variance  Reduction 

The  results  of  this  research  clearly  lead  to  the  conclusion  that  some  kind  of 
variance  reduction  techniques  applied  to  the  simulation  would  greatly  benefit  the  analyst. 
If  the  analyst  chooses  to  use  the  traditional  method  of  solving  this  kind  of  problem  (with 
only  one  realization  of  the  process)  variance  reduction  procedures  appear  to  be  critical  if 
he  hopes  to  have  any  confidence  in  tlic  solution.  If  the  analyst  chooses  to  follow  the 
approach  recommended  in  this  research  variance  reduction  will  play  a  key  role  in 
minimizing  the  number  of  extreme  points  sampled  and  aiding  in  the  comparison  between 
competing  extreme  points. 

Please  refer  to  Law  and  Kelton  (1991)  for  explanation  of  how  to  apply  variance 
reduction  techniques.  Some  techniques  that  may  be  appropriate  here  are:  multiple 
replications,  common  random  numbers,  antithetic  random  numbers,  and  control  variates. 
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9.3  Three  Step  Process 


As  stated  earlier,  the  purpose  of  this  research  was  to  investigate  the  traditional 
approach  of  solving  constrained  optimization  problems  through  simulation  and  linear 
programming.  Problems  have  been  identified  and  a  possible  solution  offered.  Figure  13 
illustrates  a  standard  estimation  of  a  response  surfaee  and  its  variance-eovariance  matrix. 


Figure  13  Step  1 


Figure  14  illustrates  a  Monte  Carlo  or  design  procedure  (using  screening)  to  sample  the 
objective  function  with  the  goal  of  including  the  tnie  optimal  extreme  point  (EP*")  in  the 
sampled  set.  It  appears  the  superior  choice  is  to  uses  the  screening  procedure  and  then 
sample  from  the  objective  function  using  a  triple  5-level  Box-Behnken  type  design. 
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Figure  14  Step  2 


Figure  15  illustrates  how  to  identify  the  "tnie"  optimal  extreme  point  (with  a  given 
probability)  and  present  the  information  to  a  decision  maker  in  both  numeric  and  visual 
form. 
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Figure  15  Step  3 


9.4  Further  Research  Areas 

This  research  could  be  continued  in  the  following  areas: 

1.  Investigate  alteniativc  sampling  designs  with  the  goal  of  decreasing  the  number  of 
design  points  and  increasing  the  "tnie"  optimal  extreme  point  sampling  accuracy. 

2.  Investigate  an  alternative  structure  to  decrease  the  number  of  design  points  in  the 
triple  modified  5-level  Box-Behnken  design. 

3.  Investigate  a  method  to  determine  a  priori  tho:  number  of  samples  that  must  be  taken  to 
achieve  a  given  confidence  level. 

4.  Investigate  an  alternative  to  Law  and  Kelton's  ranking  and  selection  procedure 
presented  in  this  research. 
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Appendix  A:  Illustration  of  Three  Step  Process 
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Appendix  C:  Box  Behnken  Type  Designs 


-1-100 
1-10  0 
-110  0 
110  0 
00-1-1 
0  0  1-1 
0  0-11 
0  0  11 
-1  0  0-1 
10  0-1 
-10  0  1 

10  0  1 

0-1-1  0 
0  1-10 
0-110 
0  110 
-10-10 
10-10 
-10  10 
10  10 
0-1  0-1 
0  10-1 
0-101 
0  10  1 

0  0  0  0 

Single  Box-Behnken  Design  (four  variables) 
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-1.0 
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0.0 

1.0 

-0.5 
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0.0 

1.0 

0.0 

0.5 

0.0 
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0.0 
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0.0 

1.0 
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Box-Behnken  Type  Design  with  5  levels  (four  variables) 
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Appendix  C;  Constants  for  the  Selection  Procedures 


p* 

'*0 

II 

II 

k  =  4 

II 

.u; 

k=:6 

II 

a 

so 

II 

k  =  9 

o 

II 

0.90 

20 

1.896 

2.342 

2.583 

2.747 

2.870 

2.969 

3.051 

3.121 

3.182 

0.90 

40 

1.852 

2.283 

2.514 

2.669 

2.758 

2.878 

2.954 

3.019 

3.076 

0.95 

20 

2.453 

2.872 

3.101 

3.258 

3.377 

3.472 

3.551 

3.619 

3.679 

0.95 

40 

2.386 

2.786 

3.003 

3.150 

3.260 

3.349 

3.422 

3.484 

3.539 

(8:606) 


Table  15.  Value  for  hi 
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(for  m  =1,  use  Table  1) 


m 

k  =  4 

k  =  5 

k=6 

k  =  7 

<x> 

II 

k  =  9 

II 

o 

P*=.90 

o 

II 

o 

2 

1.137 

1.601 

1.860 

2.039 

2.174 

2.282 

2.373 

2.450 

3 

0.782 

1.243 

1.507 

1.690 

1.830 

1.943 

2.038 

4 

0.556 

1.012 

1.276 

1.461 

1.603 

1.718 

5 

0.392 

0.843 

1.105 

1.291 

1.434 

6 

0.265 

0.711 

0.971 

1.156 

7 

0.162 

0.603 

0.861 

8 

0.075 

0.512 

9 

P*=.90 

«o=40 

N/A 

2 

1.114 

1.570 

1.825 

1.999 

2.131 

2.237 

2.324 

2.399 

3 

0.763 

1.219 

1.479 

1.660 

1.798 

1.909 

2.002 

4 

0.541 

0.991 

1.251 

1.434 

1.575 

1.688 

5 

0.381 

0.824 

1.083 

1.266 

1.408 

6 

0.257 

0.693 

0.950 

1.133 

7 

0.156 

0.587 

0.841 

8 

0.072 

0.497 

9 

II 

i/i 

no=20 

N/A 

2 

1.631 

2.071 

2.321 

2.494 

2.625 

2.731 

2.819 

2.894 

3 

1.256 

1.697 

1.952 

2.131 

2.267 

2.378 

2.470 

4 

1.021 

1.458 

1.714 

1.894 

2.033 

2.146 

5 

0.852 

1.284 

1.539 

1.720 

1.860 

6 

0.721 

1.149 

1.402 

1.583 

7 

0.615 

1.038 

1.290 

8 

0.526 

0.945 

9 

P*=.95 

No=40 

0.449 

2 

1.591 

2.023 

2.267 

2.435 

2.563 

2.665 

2.750 

2.823 

3 

1.222 

1.656 

1.907 

2.082 

2.217 

2.325 

2.415 

4 

0.990 

1.420 

1.672 

1.850 

1.987 

2.098 

5 

0.824 

1.248 

1.499 

1.678 

1.816 

6 

0.695 

1.114 

1.363 

1.541 

7 

0.591 

1.004 

1.252 

8 

0.505 

0.913 

9 

0.430 

(8:606) 


Table  16.  Value  for  h2 
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(for  m=l,  use  Table  1) 


m 

II 

k=4 

k  =  5 

k=6 

k  =  7 

k=8 

k  =  9 

0 

r-< 

II 

P*=.90 

^0=20 

2 

2.342 

2.779 

3.016 

3.177 

3.299 

3.396 

3.477 

3.546 

3 

2.583 

3.016 

3.251 

3.411 

3.532 

3.629 

3.709 

4 

2.747 

3.177 

3.411 

3.571 

3.691 

3.787 

5 

2.870 

3.299 

3.532 

3.691 

3.811 

6 

2.969 

3.396 

3.629 

3.709 

7 

3.051 

3.477 

3.709 

8 

3.121 

3.546 

9 

P*=.90 

«o=40 

3.182 

2 

2.283 

2.703 

2.928 

3.081 

3.195 

3.285 

3.360 

3.424 

3 

2.514 

2.928 

3.151 

3.302 

3.415 

3.505 

3.579 

4 

2.669 

3.081 

3.302 

3.451 

3.564 

3.653 

5 

2.785 

3.195 

3.415 

3.564 

3.675 

6 

2.878 

3.285 

3.505 

3.653 

7 

2.954 

3.360 

3.579 

8 

3.019 

3.424 

9 

P*=.95 

0 

II 

0 

3.076 

2 

2.872 

3.282 

3.507 

3.662 

3.779 

3.873 

3.952 

4.019 

3 

3.101 

3.507 

3.731 

3.885 

4.001 

4.094 

4.172 

4 

3.258 

3.662 

3.885 

4.037 

4.153 

4.246 

5 

3.377 

3.779 

4.001 

4.153 

4.269 

6 

3.472 

3.873 

4.094 

4.246 

7 

3.551 

3.952 

4.172 

8 

3.619 

4.019 

9 

P*=.95 

/io=40 

3.679 

2 

2.786 

3.175 

3.386 

3.530 

3.639 

3.725 

3.797 

3.858 

3 

3.003 

3.386 

3.595 

3.738 

3.845 

3.931 

4.002 

4 

3.150 

3.530 

3.738 

3.879 

3.986 

4.071 

5 

3.260 

3.639 

3.845 

3.986 

4.092 

6 

3.349 

3.725 

3.931 

4.071 

7 

3.422 

3.797 

4.002 

8 

3.484 

3.858 

9 

3.539 

(8:606) 


Table  17.  Value  for  [13 


55 


Appendix  D;  Main  Computer  Program  Listing 
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PROGRAM  RSMLP 


C  Written  by  ILT  R.  Garrison  Harvey  March  1,  1992 
C  This  program  investigates  the  hov?  the  noise  in  the  estimation  of 
C  a  response  surface  impacts  the  estimation  of  the  optimal  solution 
C  v7hen  the  response  surface  is  used  as  an  objective  function  of  a 
linear 
C  program. 

C 

C  This  program  is  a  research  tool,  but  V7ith  the  elimination  of 
C  unneeded  procedures  could  be  used  in  practice. 

C 

C  This  program  has  the  following  main  loops: 

C 

C  Noise  loop  (NL  times) 

C  Response  surface  (obj .  function)  sampling  (DRUMS  times) 

C  Sampling  of  the  objective  function  of  the  LP  to 

identify  the  "true"  optimal  extreme  point  (RUNS  times) 


ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

CCCCCCCCCCCC  VARIABLES  FOR  THE  LP  PORTION  CCCCCCCCCCCCCC 

INTEGER  LDA,  M,  NVAR,  RUNS,  BCDIM,  EPDIM 

PARAMETER  (M^l , NVAR=4 , LDA=M, RUNS=14  9 , BCDIM=20 , EPDIM=20 ) 

INTEGER  DP,  DRUNS,  TRUNS,  NVARY, NINT, NL, STEP 

PARAMETER  (DP=16,  DRUNS=1000,  TRUNS=RUNS*DRUNS,  NVARY=NVAR+1) 
PARAMETER  (NINT=20,  NL=5,  STEP=50) 

C  ! ! ! !  When  changing  the  parameters  remember  to  change  the  subroutines 
INTEGER  IRTYPE(M),  NOUT 

REAL  A(M,NVAR),  B(M),  C(NVAR) ,OBJTRUE(NVAR) ,  DSOL(M) 

REAL  XLB(NVAR),  XSOL(NVAR),  XUB(NVAR),  OBJ, CONSTANT 
REAL  XSOLT(NVAR),  OPTT 
EXTERNAL  DLPRS,  SSCAL,  UMACH 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
CCCC  VARIABLES  FOR  THE  MULTIVARIATE  NORMAL  GENERATION  CCC 

INTEGER  IRANK,  ISEED 

REAL  R(RUNS, NVAR) ,RSIG(NVAR, NVAR) 

EXTERNAL  CHFAC,  RNMVN,  RNSET,  UVSTA,  RNNOA,RNNOF 
EXTERNAL  CORVC,  WRRRN,  RCOV,HHSTP,  OWFRQ,  PROBP, WROPT 
EXTERNAL  SCOLR, LINRG , MRRRR 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
CCCCC  VARIABLES  FOR  THE  MONTE  CARLO  SIMULATION  CCCCCCCCC 

INTEGER  RUN, EP, BC, EPCNT(EPDIM) ,BCCNT(BCDIM) , SCREEN 
INTEGER  J,K,N,I,  SIM,  FailSamp,  SEP, S, SS, SSS,  SET 
INTEGER  CONSTRAINTS ( 3 ) 
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INTEGER  XBASIC(NL,BCDIM,M} 

REAL  NOISEMULT(NL) ,  OPTIMUM ( TRUNS, 1) 

REAL  BASIS (BCDIM, NVAR) ,  OPTBASIS(BCDIM, TRUNS) 

REAL  EXTPT(EPDIM,  NVAR),  OPTEP(EPDIM, TRUNS) 

REAL  TESTU,  TESTL,  TOL,  TEMP,  AVE,  MIN,  MAX,  S2 
REAL  OPTDR(DRUNS) 

REAL  DESIGN ( DP , NVAR+1 ) , Y ( DP, 1 ) , BHAT ( NVARY , 1 ) , VARCOV ( NVAR , NVAR) 
REAL  SEXTPT(EPDIM,NVAR) ,  box (STEP, NVAR) ,  SDEV(3) 

REAL  BA_SET(NL, BCDIM, M, NVAR) , BA (M, NVAR) , CBA(NVAR) 

REAL  CB(M),  BMAT{M,M) ,B_TEST(M) 

COMMON  OPTBASIS, OPTEP, BCCNT, BASIS, EPCNT, EXTPT, XSOL, OBJ 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

CCCCCCCCCCCCCCCC  Definitions  CCCCCCCCCCCCCCCCCCCCCCCC 

CALL  UMACH  (2,  NOUT) 

DATA  XLB/NVAR*0.0/ 

DATA  XUB/NVAR*-1.0E30/ 

DATA  A/  1 . ,  2 . ,  - 1 . , 

+  3*1., 

+  2.,  “I.,  1., 

+  l.,l.,2./ 

DATA  B/12.,  14.,  10./ 

DATA  ODJTRUE/15 . , 17 . , 18 . , 20 ./ 

CONSTANT  -  10.0 
DATA  IRTYPE/1,  1,  1/ 

C  DEFINES  THE  TYPE  OP  CONSTRAINTS 
DATA  CONSTRAINTS  /3,  0,  0/ 

C  SAMPLING  •=  1  <>  MONTE  CARLO  “2  <>  DESIGN 
SAMPLING  »  2 

C  Define  standard  deviation  multipliers 
SDEV(l)  “1.5 
SDEV(2)  “  2.75 
SDEV(3)  “4.0 

C  SCREEN  “1  <>  YES,  SCREENING,  “0  <>  NO  SCREENING 
SCREEN  =  1 

NOISEMULT(l)  =1.0 
NOISEMULT(2)  =5.0 
NOISEMOLT(3)  =9.0 
NOISEMULT(4)  =13.0 
NOISEMULT(5)  =17.0 
TOL  =  0.00001 


CCCCCCC  FULL  FACTORIAL  WITH  4  VARIABLES  CCCCCCCCCCCCCCCC 
DATA  DESIGN/16*!., 

+  "I.,!.,"!.,!.,*!.,!., "I.,!., "I.,!., "I.,!., "I.,!., “I.,!., 
+  "l.,"l.,l.,l.,“l.,"l.,l.,l.,”l.,"l.,l.,l.,”l.,”l.,l.,l.. 
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+  4*-l. ,4*1. ,4*-l. ,4*1. , 

+  8*-l. ,8*1./ 

CCCCCCC  HALF  FACTORIAL  WITH  4  VARIABLES  BLOCKING  1=1234  CCC 
C  DATA  DESIGN/8*!, 

C  +  -3 . , 1 . , 1 . , - 1 . , 1 . , - 1 . , -1 . , 1 . , 

C  +  -l.,l.,-l.,l.,-l.,l.,-l.,l., 

C  +  -1. , -1. , 1. , 1. , -1. , -1. ,1. ,1. , 

C  +  -l.,-l.,-l.,-l.,l.,l.,l.,l./ 

CCCCCCC  FULL  FACTORIAL  WITH  5  VARIABLES  CCCCCCCCCCCCCCCC 

C  DATA  DESIGN/32*!. , 

C  +  -1. , 1. , -1. ,1. , -1. , 1. , -1. ,1. , -1. , 1. , -1. , 1. , -1. ,1. , -1. ,1 

C  +  -l.,l.,-l.,l.,-l.,l.,-l.,l.,-l.,l.,-l.,l.,-l.,l.,-l.,l 

C  +  -l.,-l.,l.,l.,-l.,-l.,i.,l.,-l.,-l.,l.,l.,-l.,-l.,l.,l 

C  +  -l.,-l.,l.,l.,-l.,-l.,l.,l.,-l.,-l.,l.,l.,-l.,-l.,l.,l 

C  +  4*-l.,4*l.,4*-l.,4*l.,4*-l.,4*l.,4*-l.,4*l., 

C  +  8*-l.,8*l.,8*-l.,8*l., 

C  +16*-!., 16*1./ 

CCCCCCC  FULL  FACTORIAL  WITH  6  VARIABLES  CCCCCCCCCCCCCCCC 

C  DATA  DESIGN/64*!., 

C  +  -1. , 1. , -1. ,1. , -1. ,1. , -1. ,1., -1. ,1. , -1. ,1. , -1. ,1. , -1. ,1 

C  +  -l.,l.,-l.,l.,-l.,l.,-l.,l.,-l.,l.,-l.,l.,-l.,l.,-l.,l 

C  +  -l.,l.,-l.,l.,-l.,l.,-l.,l.,-l.,l.,-l.,l.,-l.,l.,-l.,l 

C  +  -1. ,1. , -1. ,1. , -1. , 1. , -1. ,1. , -1. ,1. , -1. , 1. , -1. , 1. , -1. , 1 

C  +  -l.,-l.,l.,l.,-l.,-l.,l.,l.,-l.,-l.,l.,l.,-l.,-l.,l.,l 

C  +  -l.,-l.,l.,l.,-l.,-l.,l.,l.,-l.,-l.,l.,l.,-l.,-l.,l.,l 

•H  -l.,-l.,l.,l.,-l.,'l.,l.,l.,-l.,-l.,l.,l.,-l.,-l.,l.,l 
+  -l.,-l.,l.,l.,-l.,-l.,l.,l.,-l.,-l.,l.,l.,-l.,-l.,l.,l 

C  +  4*-l.,4*l.,4*-l.,4*l.,4*-l.,4*l.,4*-l.,4*l., 

C  +  4*-l.,4*l.,4*-l.,4*l.,4*-l.,4*l.,4*-l.,4*l., 

C  +  8*-l.,8*l.,8*-l.,8*l.,8*-l.,8*l.,8*-l.,8*l., 

C  +  16*-1., 16*1. ,16*-!., 16*1., 

C  ^  32*-l. , 32*1./ 

C  MODIFIED  five  level  .BOX-BENKIN  DESIGN  FOR  4  VARS 


DATA  BOX/-l.,l.,-l.,l., 

4*0.0, 

+ 

-l.,l.,-l.,l.,  4*0.0,  - 

l.,l.,-l.,l.. 

5*0.0, 

+ 

-.5,  .5, -.5,  .5,  4*0.0, 

+ 

- . 5, . 5, - . 5, . 5,  4*0.0, 

.5, .5, -.5, .5, 

5*0.0, 

+ 

-.5, -.5, .5, .5, 8*0. 0,-1., 

l.,-l.,l.,4*0, 

.0, 

+ 

-l.,l.,-l.,l.,0.0. 

+ 

-l.,-l.,l.,l.,8*0.0,-.5. 

.5, -.5, .5,4*0, 

.0, 

+ 

- . 5, . 5, - . 5, .5,0.0, 

+  4*0.0,-l.,l.,-l.,l.,4*0.0,-.5,-.5, .5, .5, 

+  - .5, - .5, .5, .5, 5*0.0, 

+  4*0.0,-.5,.5,-.5, .5,4*0.0,-l.,-l.,l.,l., 

+  -l.,-l.,l.,l., 5*0.0, 
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+  4*0.0, -.5, -.5, .5, .5, -.5, -.5, .5, .5, 8*0.0, 

+  -  .  5,  ■  .  5,  .  5,  .  5, 0 . 0 , 

+  4*0. 0,-l.,-l.,l.,l.,-l.,-l.,l.,l., 8*0.0, 

+  -l.,-l.,l.,l.,0.0/ 

C  MODIFIED  three  level  BOX-BENKIN  DESIGN  FOR  5  VARS 

C  DATA  BOX/-1. ,1. , -1. ,1. ,  0.0*8,-l.,l.,-l.,l.,  0.0*12, 

C  +  -1. ,1. , -1. , 1. ,  0.0*4, -l.,l.,-l.,l.,  0.0*8, 

+  -1. , -1. ,1. , 1. ,0.0*4, -1. , 1. , -1. ,1. , 0.0*12, 

C  +  -1. , 1. , -1. , 1. , 0. 0*12, -l.,l.,-l.,l., 0.0*4, 

C 

C  +  0. 0*4, -1. , 1. , -1. , 1. , 0. 0*4, -l.,-l. 0.0*8, 

C  +  -1. , -1. , 1. , 1. , 0. 0*4, -l.,l.,-l.,l., 0.0*12, 

C 

C  +  0 .0*4, -1. , -1. , 1. , 1. , 0.0*8, -1. , -1. , 1. ,1. ,0.0*8, 

C  +  -1. , -1. ,1. , 1. ,0.0*8, -1. , -1. ,1. ,1. , 0.0*4, 

C 

C  +  0. 0*8, -1. , -1. , 1. , 1. , 0. 0*4, -l.,-l. 0.0*12, 

C  +  -l.,-l.,l.,l.,-l.,-l. 0.0*8/ 


CCCCCCCCCCCC  SOLVE  LP  FOR  TRUE  ANSWER  CCCCCCCCCCCCCCCCCC 

DO  2  K  “  1,  NVAR 
C(K)  ■=  OBJTRUE(K) 

2  CONTINUE 

C  TO  MAXIMIZE,  C  MUST  BE  MULTIPLIED  BY  -1 
CALL  SSCAL  (NVAR,  -l.OEO,  C,  1) 

C  SET  BOUNDS  ON  X  VARIABLES 
do  3  i  “  l,nvar 
xlb(i)  =  0.0 
xub(i)  -=-1.0E30 

3  continue 

CALL  DLPRS  (M, NVAR, A, M, B, B, C, IRTYPE, XLB, XUB, 

+  OPTT,XSOLT,DSOL) 

C  AFTER  RUN  DSOL  MUST  BE  MULTIPLIED  BY  -1  TO  GET  TRUE  MAX 
CALL  SSCAL  (M,  -l.OEO,  DSOL,  1) 

C  ADDING  THE  CONSTANT  TERM 

OPTT  =■  -OPTT  +  CONSTANT 

CCCCCCCCCCCCCCCCCC  END  OF  TRUE  ANSWER  FOR  LP  CCCCCCCCCCCCCCC 

CCCCCC  INITIALIZE  SEED  OF  RANDOM  NUM  GENERATOR  CCCCCCCCCCCCC 

ISEED=  1735927 
CALL  RNSET  (I SEED) 

OPEN  (UNIT  =  20,  FILE=  'final. out',  STATUS  ='NEW') 
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
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CCCCCCCCCCCCCCCC  Monte  Carlo  Simulation  CCCCCCCCCCCCCCCCCCCC 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

CCCCCCCCCCCCCCCCCCCC  NOISE  LOOP  cccccccccccccccccccccccccccc 

DO  900  N  =  1,NL 

write (*,*)  'Noise  level  ',N 

FailSamp  °  0 

EP  °  1 

BC  =■  1 

DO  4  K  =  1,EPDIM 
EPCNT(K)  “  0 
4  CONTINUE 

DO  6  K  ■=  1,BCDIM 
BCCNT(K)  -  0 
6  CONTINUE 

CCCCCCCCCCCCCCCCCC  OBJECTIVE  FUNCTION  GENERATION  LOOP  CCCCCCCCCC 
DO  500  SIM°1,DRUNS 
DO  20  J=1,DP 
Y(J,1)  =  0.0 
DO  10  K“2,NVAR+1 

Y(J,1)-  Y(J, l)+DRSIGN(J,K)*OBJTRUE(K-l) 

10  CONTINUE 

Y(J,1)'‘Y(J,1)  +NOISEMULT ( N )  * RNNOF ( )  +CONSTANT 
20  CONTINUE 

CALL  REGRESSION  (DESIGN, Y, BHAT, VARCOV) 

CCCCCCCCCCCCCC  SINGLE  SAMPLE  WITH  NO  NOISE  CCCCCCCCCCCCCCCCCCCCC 

DO  80  K  “  1,  NVAR 
C(K)  “  BHAT(K+1,1) 

80  CONTINUE 

C  TO  MAXIMIZE,  C  MUST  BE  MULTIPLIED  BY  -1 
CALL  SSCAL  (NVAR,  -l.OEO,  C,  1) 

C  Solve  LP 

CALL  DLPRS  (M, NVAR, A, M, B, B, C, IRTYPE, XLB, XUB, 

+  OBJ,XSOL,DSOL) 

C  AFTER  RUN  DSOL  MUST  BE  MULTIPLIED  BY  - 1  TO  GET  TRUE  MAX 
CALL  SSCAL  (M,  -l.OEO,  DSOL,  1) 

C  ADDING  THE  CONSTANT  TERM 
OBJ  =  -OBJ  +  Bhat(l,l) 


CCCCCCCCCCCC  END  SINGLE  SAMPLE  WITH  NO  NOISE  CCCCCCCCCCCCCCCCCCC 
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C  OBTAIN  THE  CHOLESKY  FACTORIZATION  (GEN  OF  MULTIVARIATE  DIST) 

CALL  CHFAC  (NVAR, VARCOV, NVAR, 0 . 00001, IRANK, RSIG, NVAR) 

CALL  RNMVN  (RUNS,  NVAR,  RSIG,  NVAR,  R,  RUNS) 

C  RNMVN  (NR,  K,  RSIG,  LDRSIG,  R,  LDR) 

C  NR  -  if  RANDOM  MULTIVARIATE  NORMAL  VECTORS  TO  GENERATE  (INPUT) 

C  K  -  (EQ  NVAR)  LENGTH  OF  THE  MULTIVARIATE  NORMAL  VECTOR  (INPUT) 

C  RSIG  -  UPPER  TRIANG  MATRIX,  K  BY  K,  CONTAINING  THE  CHOLESKY  FACTOR 
FOR  THE 

C  VARIANCE- COVARIANCE  MATRIX  (INPUT) 

C  LDRSIG  -(RUNS)  LEADING  DIM  OF  RSIG  EXACTLY  AS  SPECIFIED  IN  THE 
CALLING 
C  PROGRAM 

C  (INPUT) 

C  R  -  NR  BY  K  MATRIX  CONTAINING  THE  RANDOM  VECTOR  IN  ITS  ROWS 

C  LDR  (RUNS)  -  LEADING  DIM  OF  R  EXACTLY  AS  SPECIFIED  IN  THE  DIMENSION 

C  STATEMENT  OF  THE  CALLING  PROGRAM  (INPUT) 

CCCCCCCCCCCCCCCC  OBJECTIVE  FUNCTION  SAMPLING  LOOP  CCCCCCCCCCC 

OPTDR(SIM)“0 .0 
SEP  ■=•  1 
SSS  “0 
SET  “  0 

DO  390  RUN  =  1,  RUNS 

DO  100  K  *■  1,  NVAR 
IF  ( SAMPLING. EQ.l)  THEN 
C(K)  =  DHAT(K+1,1)  +  R(RUN,K) 

ELSE 

SS  =  RUN  -  SSS*STEP 
S  “  SSS  +1 
SSS“INT( RUN/STEP) 

C(K)“BHAT(K+1, 1)+ 

+  BOX(SS,K)*SDEV(S)*SQRT(VARCOV(K,K) ) 


END  IF 

100  CONTINUE 

IF  (SCREEN. EQ.l)  GOTO  120 

C  TO  MAXIMIZE,  C  MUST  BE  MULTIPLIED  BY  -1 
CALL  SSCAL  (NVAR,  -l.OEO,  C,  1) 

C  Solve  LP 

CALL  DLPRS  (M, NVAR, A, M, B, B, C, IRTYPE, XLB, XUB, 

+  OBJ,XSOL,DSOL) 

C  AFTER  RUN  DSOL  MUST  BE  MULTIPLIED  BY  -1  TO  GET  TRUE  MAX 
CALL  SSCAL  (M,  -l.OEO,  DSOL,  1) 

C  ADDING  THE  CONSTANT  TERM 

OBJ  =  -OBJ  +  Bhat(l,l) 
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C  DLPRS  (M,  NVAR,  A,  LDA,  BL,  BU,  C,  IRTYPE,  XLB,  XUB,  OBJ,  XSOL, 
DSOL) 

C  M  -  i  OF  CONSTRAINTS  (INPUT) 

C  NVAR  -  S  OF  VARIABLES  (INPUT) 

C  A  -  MATRIX  OF  DIM  m  BY  NVAR  CONTAINING  THE  COEFFICIENTS  OF  THE  M 
CONST . 

C  (INPUT) 

C  LDA  -  LEADING  DIM  OF  A  EXACTLY  AS  SPECIFIED  IN  THE  DIM  STATEMENT 
(INPUT) 

C  BL,  BU  -  UPPER  S  LOWER  BOUNDS  VECTOR  ENGTH  M  (INPUT) 

C  C  -  VECTOR  LENGTH  NVAR  •=■  COEFF.  OF  OBJ  FUNCT  (INPUT) 

C  IRTYPE  -  VECTOR  LENGTH  M  °  TYPE  OF  CONSTRAINTS  EXCLUSIVE  OF  SIMPLE 
BOUNDS 

C  WHERE  IRTYPE(I)“  0,1, 2, 3  INDICATES  .EQ.,  .LE.,  .GE.,  AND  RANGE 
CONSTRAINTS 
C  RESPECTIVILY 

C  XLB,XUB  -  VECTORS  LENGTH  NVAR  LOWER  &  UPPER  BOUNDS  OF  VARIABLES 
C  IF  NO  BOUNDS  THEN  SET  XLB  1.0E30,  OR  XUB  “  -1.0E30 
C  DEPENDING  WHICH  IS  UNBND 

C  OBJ  -  VALUE  OF  OBJECTIVE  FUNCTION  (OUTPUT) 

C  XSOL  -  VECTOR  LENGTH  NVAR  -  PRIMAL  SLOUTION  (OUPUT) 

C  DSOL  -  VECTOR  LENGTH  M  “  DUAL  SOLUTION  (OUTPUT) 

GOTO  155 

C  selection  procedure  to  follow 

120  IP  (RUN.GT.l)  THEN 

DO  135  J»  1,SET 

C  DEFINE  BASIC  COEFFICIENTS 
DO  121  1=1, M 
TEMP  =  XBASIC(N, J, I) 

IF  (TEMP.EQ.O)  THEN 
CB(I)  “  0.0 
ELSE 

CB(I)  =  C(XBASIC(N, J, I) ) 

ENDIF 

121  CONTINUE 

C  CALCULATE  Cb*Binv*A  VECTOR  GIVEN  BA  =  Binv*A 
DO  124  K=  1,  NVAR 
CBA(K)  =  0.0 
DO  122  I  =  1,  M 

CBA(K)=CBA(K)+CB(I)*BA_SET(N, J,I,K) 

122  CONTINUE 

124  CONTINUE 

C  OPTIMALITY  TEST  -  BRANCH  IF  NOT  OPTIMAL 
DO  130  1=  1,NVAR 

IF(C(I) -CBA(I) .GT.tol)  GOTO  134 
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130 


CONTINUE 


C  BASIS  ALREADY  SAMPLED  -  OPTIMAL  CONDITIONS  MET,  DON'T  SAMPLE 
C  Branch  and  sample  objective  function  agian 
GOTO  390 

134  temp=0 

135  CONTINUE 

ELSE 

SET  =  0 
ENDIF 

SET  =■  SET  +1 

C  TO  MAXIMIZE,  C  MUST  BE  MULTIPLIED  BY  -1 
CALL  SSCAL  (NVAR,  -l.OEO,  C,  1) 

C  Solve  LP 

CALL  DLPRS  (M, NVAR, A, M, B, B, C, IRTYPE, XLB, XUB, 

+  OBJ,XSOL,DSOL) 

C  AFTER  RUN  DSOL  MUST  BE  MULTIPLIED  BY  -1  TO  GET  TRUE  MAX 
CALL  SSCAL  (M,  -l.OEO,  DSOL,  1) 

C  ADDING  THE  CONSTANT  TERM 

OBJ  “  -OBJ  +  Bhat(l,l) 

COUNT  »  0 

DO  137  I  -  1,  NVAR 

IF{XSOL(l) .NE.O.O)  THEN 
COUNT  -  COUNT+1 
XBASIC(N, SET, COUNT)  =  I 
DO  136  J  ■=>  1,  M 

BMAT(J, COUNT)  =  A(J,I) 

136  CONTINUE 
ENDIF 

137  CONTINUE 

IF(COUNT.GT.M)  WRITE  (*,*)  'ERROR  COUNT  >  M' 

C 

CCC  FOLLWING  IS  EXECUTED  WHEN  A  BASIC  VARIABLE  IS  NOT  A  DECISION 
VARIABLE  CCC 
C 

IF  (COUNT. LT.M)  THEN 
DO  142  I  =  COUNT+1,  M 
XBASIC(N,SET,I)=  0.0 
DO  140  J  =  1,  M 
BMAT(J,I)  =0.0 
140  CONTINUE 

142  CONTINUE 

DO  148  1=  1,  M 
B_TEST(I)  =0.0 
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DO  146  J  =  1,  NVAR 

B_TEST ( I ) =B_TEST ( I ) +XSOL( J ) * A ( I , J ) 

146  CONTINUE 

148  CONTINUE 

DO  150  I  =■  1,  M 

IF  ( (B_TEST(I)-B(I) .NE.0.0) .AND. ({DSOL(I) .LE.TOL) 

+  .AND. (DSOL(I) .GE. -TOL) ) )  THEN 
COUNT=*COUNT+l 
BMAT(I,COUNT)  »  1 
ENDIF 

150  CONTINUE 
ENDIF 

IF  (COUNT. GT.M)  WRITE  (*,*)  'ERROR;  TOO  MANY  B  coin  DEFINED' 

C  INVERTING  B 

CALL  LINRG(M,BMAT,M,BMAT,M) 

C  FINDING  Binv*A 

CALL  MRRRR ( M , M , BMAT , M , M , NVAR , A , M, M, NVAR , BA , M ) 

CCCCCCCC  END  OF  FINDING  B  WHEN  BASIC  VAR  IS  NOT  A  DECISION  VAR 
CCCCCCCCCCCC 

DO  153  J-1,M 

DO  151  I»  1,NVAR 

BA_SET(N,SET, J,I)  “  BA(J,1) 

151  CONTINUE 

153  CONTINUE 

C  DEFINE  RUNNING  STATISTICS 

155  Ol’TDR(SIM)  =  OPTDR(SIM)  +  OBJ 

OPTIMUM(RUN+  RUNS* (SIM-1), 1)  “  OBJ 

CCCCCCCC  DEFINE  FOR  PER  SAMPLE  TESTING  CCCCCCCCCCCCCCCCCCCC 
IF  (RUN.EQ.l)  THEN 
DO  157  K  °  1,  NVAR 
SEXTPT(1,K)  =  XSOL(K) 

157  CONTINUE 

ENDIF 


C  Define  first  Extreme  pt  &  Decision  Set  for  each  run 

IF  ( (RUN.EQ.l) .AND. (SIM. EQ.l))  THEN 
BCCNT(l)  =  1 
OPTBASIS(l,l)  =  OBJ 
OPTEP(l,l)  =  OBJ 
EPCNT(l)  =  1 
DO  158  K  =  1,  NVAR 
EXTPT(1,K)  =  XSOL(K) 

BASIS(1,K)  =  XSOL(K) 

158  CONTINUE 
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GO  TO  390 
ENDIF 

CCCCCCCCCCCCC  TEST  COMPARE  EXTREME  POINTS  PER  SAMPLE  CCCCCC 
J  =  SEP 

160  IF  (J.GT.O)  THEN 

DO  170  K  =  1,  NVAR 
TESTL  =  XSOL(K)  -  TOL 
TESTU  -  XSOL(K)  +  TOL 

IF  ( (SEXTPT(J,K) .LT. TESTL) .OR. (SEXTPT(J,K) .GT. TESTU) )  THEN 
J“  J-1 

IF  (J.EQ.O)  THEN 
GO  TO  175 
ELSE 

GO  TO  160 
ENDIF 
ENDIF 

170  CONTINUE 

GO  TO  199 

C  DEFINE  A  NEW  EXTREME  POINT  FOR  THIS  RUN  ONLY 

175  SEP  “  SEP  +1 

DO  180  K  -  1,NVAR 

SEXTPT(SEP,K)  “  XSOL(K) 

180  CONTINUE 

ENDIF 

CCCCCCCCCCCCCC  TEST  TO  COMPARE  EXTREME  POINTS  CCCCCCCCCCCC 

199  J  ”  EP 

200  IF  (J.GT.O)  THEN 

DO  250  K  =  1,  NVAR 

TESTL  =  XSOL(K)  -  TOL 

TESTU  “  XSOL(K)  +  TOL 

IF  ( (EXTPT(J,K) .LT. TESTL) .OR. (EXTPT(J,K) .GT. TESTU) )  THEN 
J“  J-1 

IF  (J.EQ.O)  THEN 
GO  TO  275 
ELSE 

GO  TO  200 
ENDIF 
ENDIF 

250  CONTINUE 

CALL  ASSIGNEP(J) 

J  =  J  -1 
GO  TO  290 

C  DEFINE  A  NEW  EXTREME  POINT 

275  EP  =  EP  +1 

CALL  ASSIGNEP(EP) 

ENDIF 
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CCCCCCCC  TEST  TO  COMPARE  DECIOSON  VAR  SET  CHANGES  CCCCCCCC 
290  J  =  BC 

300  DO  350  K  “  1,  NVAR 

IF  ( ( ( (BASIS(J,K) .LE.TOL) .AND. (BASIS(J,K) .GE. -TOL) ) 

+  .AND. ( (XSOL(K) .GE.TOL} .  OR. (XSOL(K) . LE. -TOL) ) ) 

+  .OR. (( (BASIS(J,K) .GE.TOL) .  OR .( BASIS (J, K) . LE . -TOL) ) ) 

+  .AND. ( (XSOL(K) .LE.TOL) .AND. (XSOL(K) .GE. -TOL) ) )  THEN 
GO  TO  355 
ENDIF 

350  CONTINUE 

CALL  ASSIGNBC  (J) 

GO  TO  390 

355  IP  (J.LE.l)  THEN 

BC  “  BC+1 
CALL  ASSIGNBC  (BC) 

ELSE 

J  “  J  -1 
GO  TO  300 
ENDIF 

390  CONTINUE 

CCCCCCCC  TEST  IF  TRUE  EXTREME  POINT  WAS  SAMPLED  CCCCCCCCCCC 
J  “  SEP 

400  IF  (J.GT.O)  THEN 

DO  470  K  -  1/  NVAR 
TESTL  “  XSOLT(K)  -  TOL 
TESTU  •»  XSOLT(K)  +  TOL 

IF  ( (SEXTPT(J,K) .LT. TESTL) .OR. (SEXTPT(J,K) .GT. TESTU ) )  THEN 
J  “  J-1 

IF  (J.EQ.O)  THEN 
GO  TO  480 
ELSE 

GO  TO  400 
ENDIF 
ENDIF 

470  CONTINUE 

GO  TO  500 

C  DEFINE  A  NEW  EXTREME  POINT  FOR  THIS  RUN  ONLY 

480  FailSamp  =  FailSamp  +  1 
ENDIF 

500  CONTINUE 


cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

CCCCCCCCCCCCCC  END  IF  FOR  MONTE  CARLO  SIM  CCCCCCCCCCCCCCCCC 
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 


cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

CCCCCCCCCCCC  CALCULATE  RESULTS  AND  PRINT  OUTPUT  CCCCCCCCCCCC 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

CALL  UMACH(-2,20) 

IF  ((XLB(l).NE.  0. ) -OR. (XLB(2) .NE.  0.))  THEN 

WRITE(20, *)' *******  WARNING  WARNING  ************' 

WRITE  (20,*)  '  LOWER  BOUNDS  ON  X  CORRUPT' 

WRITE  (20,*)  'XLB  =  ' , XLB 
ENDIF 

WRITE  (20,*)  'NOISE  MULTIPLIER  (SD)  FOR  NORMAL  NOISE:  ', 

+  NOISEMULT(N) 

TEMP  =  NOISEMULT(N)*(l/real(DP))**.5 

WRITE  (20,*)  'STANDARD  ERROR  IS  :',TEMP 

if  (screen. eq. 1)  WRITE(20,*)  'This  is  a  screening  run' 

if  ( sampling. eq. 2)  write (20,*)  '  Using  a  Design  to  sample' 

WRITE  (20,*)  '  jj  Times  failed  to  sample  true  extreme  pt', 

+  FailSamp 

write (20,*)  'Standard  deviation  mult  set  =',SDEV 
IF  ( FailSamp. GT. 0)  then 

TEMP  °  100*REAL(FailSamp)/real(DRUNS) 

WRITE  (20, '(A,  F7.3)')  '%  failures  overall:  ',  TEMP 
ENDIF 

WRITE  (20,*) 'Number  of  Objective  Function  Samples:  ',DRUNS 
WRITE  (20,*) 'Number  of  Runs  per  Obj  Function:  ',  RUNS 
WRITE  (20,*) 'Total  Number  of  Points  Tested:  ',TRUNS 
WRITE  (20,*) 'The  True  Objective  Function:  ' 

WRITE  (20,*)  'const' ,CONSTANT, '+' ,OBJTRUE 

WRITE  (20,*)  'Sample  Generated  Objective  Function' 

WRITE  (20,*)  BHAT 

CALL  WRRRN  ('Constraint  Matrix ', LDA, NVAR, A, LDA, 0 ) 

WRITE  (20, *)  '  ' 

WRITE  (20,*)  'The  RHS  is:  ',  B 

WRITE  (20,*)  ' *True  Optimal  Answer:  ',  OPTT 

WRITE  (20,*)  ' *True  optimal  Extreme  Point:' 

WRITE  (20,*)  XSOLT 

CALL  WRRRN  ('Design  Matrix ', DP, NVAR+1, DESIGN, DP, 0 ) 

WRITE  (20,*)  'Sample  response  vai'aible  Y:' 

WRITE  (20,*)  Y 

CALL  WRRRN  ('Sample  Variance-Covariance  Matrix', 

+  NVAR, NVAR, VARCOV, NVAR, 0) 

CALL  WRRRN  ('Sample  Cholesky  Factorization  Matrix', 

+  NVAR, NVAR, RSIG, NVAR, 1) 
if  (sampling .eq. 2)  then 
CALL  WRRRN  ('SAMPLING  DESING  (BOX-BEHNKEN) ' , 

+  step, NVAR, BOX, step, 0) 
end  if 

WRITE  (20, *)  ' ' 

WRITE  (20,*) 

+  '!!!!!!!!!!!!!!!!!!  DECISION  VARIABLES  !!!!!!!!!!!!!!!!!!' 

WRITE  (20,*)  'S  OF  DECISION  VARIABLE  SET  CHANGES  :  ' , BC 
WRITE  (20, *)  ' ' 
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605 


608 


620 


680 


DO  680  K  =  1,  BC 
WRITE  (20,*)  ' 

WRITE  (20,*)  'Decision  var  set  S  ',K 

DO  605  J  =  1,  NVAR 

WRITE  (20, ' (2X,fl0.4) ')  BASIS(K,J) 

CONTINUE 

WRITE  (20,*)  'if  of  occurances  of  this  basis:  ',BCCNT(K) 
if  (sampling .eq. 1)  then 

TEMP  =  (100*  REAL(BCCNT(K) ) )/TRUNS 
else 

TEMP  =  (100*REAL(BCCNT(K) ) )/DRUNS 
ENDIF 

WRITE  (20, '(A,  F7.3)')  '%  of  overall  occurance  is:  ',  TEMP 
AVE  =0.0 

DO  608  J  =  1,  BCCNT(K) 

AVE  =  AVE  +  OPTBASIS(K, J) 

CONTINUE 

AVE  =  AVE/REAL(BCCNT(K) ) 

WRITE  (20, '(A,  F10.4)')  'AVE.  OPTIMUM  IS:  ',AVE 

WRITE  (20,  '  (A,F10.4)  '  ) 'Bias  opt  est  (Ave  Opt-Tr.ue  Opt):  ', 

+  AVE  -  OPTT 
S2  =  0.0 

MIN  =  OPTBASIS(K, 1) 

MAX  =  MIN 

DO  620  J  =  1,  BCCNT(K) 

S2  =  S2+  (AVE  -  OPTBASIS(K, J) )**2 

IF  (MAX .LT.OPTBASIS(K,J) )  THEN 
MAX  ■=  OPTBASIS(K, J) 

ENDIF 

IF  (MIN.GT.OPTBASIS(K,J) )  THEN 
MIN  =  OPTBASIS  (K,J) 

ENDIF 

CONTINUE 

S2  =  S2/BCCNT(K) 

WRITE  (20, ' (A,  F13,4)')  'The  population  variance  is:  ',  S2 
WRITE  (20, ' (A,  F13.4)')  'The  maximum  value  is:  ',  max 
WRITE  (20, '(A,  F13.4)')  'The  minimum  value  is:  ',  min 

CONTINUE 

WRITE  (20,*)  ' ' 

WRITE  (20,*)  '  ' 

WRITE  (20,*) 
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+  '!!!!!!!!!!!!!!!!!!!  THE  EXTREME  POINT  !!!!!!!!!!!!!!!!!!!' 
WRITE  (20,*)  'NUMBER  OF  EXT  POINTS  VISITED  IS:  EP 

DO  770  K  =  1,  EP 
WRITE  (20,*)  ' 

WRITE  (20,*)  'Extreme  Point  g  ',  K 
DO  715  J  =  1,  NVAR 

WRITE  (20, ' (2X,fl0.4) ' )  EXTPT(K,  J) 

715  CONTINUE 

WRITE  (20,*)  'g  OF  EXT  PT  VISITS  ARE;  ',EPCNT(K) 
if  (sampling .eq. 1)  then 

TEMP  =  (100*  REAL(EPCNT(K) ) )/TRUNS 
else 

TEMP  =  (100*REAL(EPCNT(K) ) )/DRUNS 
ENDIF 

WRITE  (20, '(A,  F7.3)')  '%  OF  OVERALL  IS:  ',  TEMP 
AVE  =0.0 

DO  718  J  =  1,  EPCNT(K) 

AVE  =■  AVE  +  OPTEP(K,J) 

718  CONTINUE 

AVE  =  AVE/REAL(EPCNT(K) ) 

WRITE  (20, '(A,  F10.4)')  'AVE.  OPTIMUM  IS :  ',AVE 
WRITE  (20, '(A,  F10.4) ' ) 'Bias  opt  est  (Ave  Opt-True  Opt):  ', 
+  AVE  -  OPTT 

C  CALCULATE  TRUE  Z*  BASED  ON  THIS  EXTREME  POINT 
WRITE(20,*)  '' 

TEMP  =0.0 
DO  719  J=1,NVAR 

TEMP  =  TEMP  +  OBJTRUE(J)*BASIS(K, J) 

719  CONTINUE 

TEMP  =  TEMP  +  CONSTANT 

WRITE  (20, '(A,  F10.4)')  'TRUE  Z*  WITH  TRUE  C  IS:  ' , TEMP 
WRITE  (20, '(A,  F10.4)')  'TRUE  BIAS  (Z*  -  Z  optimal):  ', 

+  TEMP  -  OPTT 

WRITE  (20, ' (A,  F10.4) ' ) 

4  'Difference  between  expected  optimal  and  true' , TEMP -AVE 


C  CALCULATE  VARIANCE,  MINIMUM  &  MAXIMUM 
S2  =  0.0 

MIN  =  OPTEP(K,l) 

MAX  =  MIN 

DO  741  J  =  1,  EPCNT(K) 

S2  =  S2+  (AVE  -  OPTEP(K, J) )**2 

IF  (MAX.LT.OPTEP(K, J) )  THEN 
MAX  =  OPTEP(K,J) 

ENDIF 

IF  (MIN.GT.OPTEP(K, J) )  THEN 
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MIN  =  OPTEP  (K,J) 

ENDIF 

741  CONTINUE 

S2  =  S2/EPCNT(K) 

V7RITE  (20, '(A,  F13.4)')  'The  population  variance  is:  S2 

WRITE  (20, '(A,  F13.4)’)  'The  maximum  value  is:  ',  max 
WRITE  (20, '(A,  F13.4)')  'The  minimum  value  is:  ',  min 
770  CONTINUE 

WRITE  (20,*)  " 

WRITE  (20,*)  '' 

WRITE  (20,*) 

+  '!!!!!!!!!!!!!!!!!!!!  OVERALL  RESULTS  !!!!!!!!!!!!!!!!!!!!!!' 
WRITE  (20,*)  '' 

AVE  =0.0 

DO  800  J  =  1,  TRUNS 

AVE  =  AVE  +  OPTIMUM (J,l) 

800  CONTINUE 

AVE  =  AVE/TRUNS 

WRITE  (20, '(A,  F10.4)')  'Overall  Mean  Optimum  is:  ',AVE 
WRITE  (20, '(A,  FIO . 4 )')' Overall  Bias  (Ave  Opt-True  Opt):  ', 

+  AVE  -  OPTT 
S2  =  0.0 

MIN  =  OPTIMUM (1,1) 

MAX  =  MIN 

DO  820  J  =  1,  TRUNS 

S2  =  S2+  (AVE  -  OPTIMUM(J, 1) )**2 

IF  (MAX.LT.OPTIMUM(J,l) )  THEN 
MAX  =  OPTIMUM (J,l) 

ENDIF 

IF  (MIN.GT. OPTIMUM ( J, 1) )  THEN 
MIN  =  OPTIMUM(J,l) 

ENDIF 

820  CONTINUE 

S2  =  S2/ (TRUNS- 1) 

WRITE  (20, '(A,  F13.4)')  'The  overall  sample  variance  is:  ',  S2 
WRITE  (20, '(A,  F13.4)')  'The  overall  maximum  value  is:  ',  max 
WRITE  (20, '(A,  F13.4)')  'The  overall  minimum  value  is:  ',  min 
WRITE  (20,*)  ' ' 

C 

WRITE  (20,*) ' ' 

WRITE (20,*) 'AVE  OPTIMAL  PER  OBJECTIVE  FUNCTION' 

DO  855  J  =  1,  DRUNS 

OPTDR(J)  =  OPTDR(J)/RUNS 
855  CONTINUE 
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AVB  =0.0 

DO  856  J  =  1,  DRUMS 

AVE  =  AVE  +  OPTDR(J) 

856  CONTINUE 

AVE  =  AVE/DRUNS 

WRITE( 20 , ' ( A,  FIO . 4 ) ' ) 'Mean  of  Mean  Opt  per  obj  funct:  ',AVE 

S2  =  0.0 
MIN  »  OPTDR(l) 

MAX  “  MIN 

DO  85?  J  »  1,  DRUMS 

S2  =  S2+  (AVE  -  OPTDR(J))**2 

IF  (MAX.LT.OPTDR(J) )  THEN 
MAX  -  OPTDR{J) 

ENDIF 

IF  (VIN.GT.OPTDR(J) )  THEN 
MIN  =  OPTDR(J) 

ENDIF 

857  CONTINUE 

S2  -  S2/(DRUNS-1) 

WRITE  (20, '{A,  F13.4)')  'The  sample  var(mean  opt):  S2 

WRITE  (20, '(A,  F13.4)')  'The  maximum  value  is:  ',  max 
WRITE  (20, '(A,  F13.4)')  'The  minimum  value  is:  ',  min 
WRITE  (20,*)  ' ' 

DO  860  J  -  1,15 
WRITE  (20,*)  '■ 

860  CONTINUE 

900  CONTINUE 

CLOSE  (20) 


END 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
CCCCCCCCCCCCCCCCCCCC  END  t’ROGRAM  CCCCCCCCCCCCCCCCCCCCCCCCCC 
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

SUBROUTINE  ASSIGNBC  (TEMP) 

INTEGER  TEMP,  KK,  NVAR,  BCDIM,  EPDIM,RUNS 
PARAMETER  (NVAR=4  ,  BCDIM=20,  EPDIM=20,  RUNf. --149  ) 

INTEGER  DP,  DRUMS,  TRUNS,  NVARY 

PARAMETER  (DP“16,DRUNS=1000,TRUNS=RUNS*DRUNS, NVARY=NVAR+1) 

INTEGER  EPCNT(EPDIM) ,  BCCNT(DCDIM) 

REAL  OPTBASIS' BCDIM, TRUNS) 

REAL  OPTEP(EPDIM, TRUNS) ,BASIS(BCDIM, NVAR) 
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REAL  XSOL(NVAR) ,OBJ,  EXTPT(EPDIM, NVAR) 

COMMON  OPTBASIS, OPTEP, BCCNT, BASIS, EPCNT, EXTPT, XSOL, OBJ 

C  COUNT  OCCURANCES  OF  EACH  DECISION  SET  CHANGE 
BCCNT (TEMP)  «  BCCNT (TEMP)  +  1 

C  DEFINE  DECISION  SET 

OPTBASIS (TEMP, BCCNT ( TEMP ) )  »  OBJ 
DO  1010  KK  “  1,  NVAR 

BASIS (TEMP, KK)  “  XSOL(KK) 

1010  CONTINUE 

RETURN 

END 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

SUBROUTINE  ASSIGNED  (TEMP) 

INTEGER  TEMP,  KK,  NVAR,  BCDIM,  EPDIM, RUNS 
PARAMETER  ( NVAR=>4  ,  BCDIM°20 ,  EPDIM“20 ,  RUNS=i4  9  ) 

INTEGER  DP,  DRUNS,  TRUNS,  NVARY 

PARAMETER  ( DP=«16 , DRUNS-IOOO , TRUNS-RUNS*DRUNS , NVARY-NVAR+1) 

INTEGER  EPCNT(EPDIM) ,  BCCNT(BCDIM) 

REAL  OPTBASIS (BCDIM, TRUNS) 

REAL  OPTEP ( EPDIM, TRUNS ) , BASIS ( BCDIM, NVAR ) 

REAL  XSOL (NVAR), OBJ,  EXTPT(EPDIM, NVAR) 

COMMON  OPTBASIS, OPTEP, BCCNT, BASIS, EPCNT, EXTPT, XSOL, OBJ 

C  COUNT  OCCURANCES  OF  EACH  EXTREME  POINT 
EPCNT (TEMP)  “  EPCNT (TEMP)  +  1 

C  DEFINE  EXTREME  POINT 

OPTEP(TEMP, EPCNT ( TEMP ) )  =  OBJ 
DO  1020  KK  =  1,  NVAR 

EXTPT ( TEMP, KK)  =  XSOL(KK) 

1020  CONTINUE 

RETURN 

END 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

SUBROUTINE  REGRESS ION ( X, Y, BHAT, VARCOV) 

INTEGER  NVAR,  NVARY,DP,KK,  JJ 
PARAMETER  ( NVAR“4 , NVARY=NVAR+1, DP=16 ) 

REAL  X(DP,NVARy) ,Y(DP, 1 ) , BHAT ( NVARY , 1) , VARCOV(NVAR, NVAR) 
REAL  XX ( NVARY , NVARY ) , C ( NVARY , 1 ) , INV ( NVARY , NVARY ) 

REAL  MSE,  EY(DP) 

r  ETERNAL  ^KTXF,LSGRR,MXTYF,MRRRR 
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C  COMPUTE  (X'X) 

CALL  MXTXF ( DP , NVARY , X , DP , NVARY , XX , NVARY ) 

C  INVERT  THE  MATRIX 

CALL  LSGRR ( NVARY , NVARY , XX , NVARY ,  0 . 0  0  0  0 1 , 2 , 1 NV , NVARY ) 

C  COMPUTE  X'Y 

CALL  MXTYF ( DP , NVARY , X , DP , DP , 1 , Y , DP , NVARY , 1 , C , NVARY ) 

CALL  MRRRR( NVARY, NVARY, INV, NVARY, NVARY, 1, 

+  C , NVARY , NVARY , 1 , BHAT , NVARY ) 

C  CALCULATE  MSE 
MSE  "0.0 

DO  1105  KK  “  1,DP 
EY{KK)  “  0.0 
DO  1.104  JJ  “  1,  NVARY 

EY(KK)“EV(KK)+X(KK,  JJ)*BI!AT(JJ,  1) 

1104  CONTINUE 

1105  CONTINUE 

DO  1107  KK“  1,  DP 

MSE  *■  MSE  I-  (Y(KK,1)-EY(KK))**2 
1107  CONTINUE 

MSE  -  MSE/(  DP  -  NVARY) 

C  CALCULATE  VARIANCE- COVARIANCE  MATRIX 
DO  1110  K“2, NVARY 

DO  1100  J“2, NVARY 

VARCOV (J-1,K-1)"INV(J,K)*MSE 
1100  CONTINUE 

1110  CONTINUE 

RETURN 

END 

cccccccccccccccccccccccccccccccccccccccccccccccccccc 
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Appendix  E:  Sample  Computer  Program  Output 
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NOISE  MULTIPLIER  (SD)  FOR  NORMAL  NOISE:  1.00000 

STANDARD  ERROR  IS  :  0.250000 

This  is  a  screening  run 

Using  a  Design  to  sample 
H  Times  failed  to  sample  true  extreme  pt  0 
Standard  deviation  mult  set  =  1.50000  2.75000 

4.00000 


Number 

of  Objective  Function 

Samples : 

1000 

Number 

of  Runs  per  Obj  Function:  149 

Total  Number  of  Points  Tested 

:  149000 

The  True  Objective  Function: 

const 

10.00000+  15.0000 

17.0000  18. 

,0000 

20.0000 

Sample 

Generated  Objective  Function 

9.5 

3095  15.0707  17. 

1053 

17.7075 

19.9385 

Constraint  Matrix 

1  2 

3 

4 

1 

1  1 

2 

1 

2 

2  1 

- 

1 

1 

3 

-1  1 

1 

2 

The  RHS 

is:  12.0000  14 

.0000 

10.00000 

*True  Optimal  Answer:  216 

.000 

*True  optimal  Extreme  Point: 

2.00000  8.00000  0. 

2.00000 

Design 

Matrix 

1  2 

3 

4 

5 

1 

1  -1 

-1 

-1 

-1 

2 

1  1 

-1 

-1 

-1 

3 

1  -1 

1 

-1 

-1 

4 

1  1 

1 

-1 

-1 

5 

1  -1 

-1 

1 

-1 

6 

1  1 

-1 

1 

-1 

7 

1  -3 

1 

1 

-1 

8 

1  1 

1 

1 

-1 

9 

1  -1 

-1 

-1 

1 

10 

1  1 

-1 

-1 

1 

11 

1  -1 

1 

-1 

1 

12 

1  1 

1 

-1 

1 

13 

1  -1 

-1 

1 

1 

14 

1  1 

-1 

1 

1 

15 

1  -1 

1 

1 

1 

16 

1  1 

1 

1 

1 

Sample 

response  varaible  Y: 

-59. 

8876  -28.8743  -25. 

9255 

4.72923 

-25.1038 

4.60987 

9.42563 

37. 

7662  -20.6727  9.03708 

12.6977 

43.4840 

14.8402 

i  45.4566 

50.3080  80.6046 


Sample  Variance -Covariance  Matrix 

12  3  4 
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1 

0.05979 

0.00000 

0.00000 

0.00000 

2 

0.00000 

0.05979 

0.00000 

0.00000 

3 

0.00000 

0.00000 

0.05979 

0.00000 

4 

0.00000 

0.00000 

0.00000 

0.05979 

Sample  Cholesky  Factorization  Matrix 

1 

2 

3 

4 

1 

0.2445 

0.0000 

0.0000 

0.0000 

2 

0.2445 

0.0000 

0.0000 

3 

0.2445 

0.0000 

4 

0.2445 

SAMPLING  DESING 

(BOX-BEHNKEN) 

1 

2 

3 

4 

1 

-1.0 

-0.5 

0.0 

0.0 

2 

1.0 

-0.5 

0.0 

0.0 

3 

-1.0 

0.5 

0.0 

0.0 

4 

1.0 

0.5 

0.0 

0.0 

5 

0.0 

0.0 

-1.0 

-0.5 

6 

0.0 

0.0 

1.0 

-0.5 

7 

0.0 

0.0 

-1.0 

0.5 

8 

0.0 

0.0 

1.0 

0.5 

9 

-1.0 

0.0 

0.0 

-0.5 

10 

1.0 

0.0 

0.0 

-0.5 

11 

-1.0 

0.0 

0.0 

0.5 

12 

1.0 

0.0 

0.0 

0.5 

13 

0.0 

-1.0 

-0.5 

0.0 

14 

0.0 

1.0 

-0.5 

0.0 

15 

0.0 

-1.0 

0.5 

0.0 

16 

0.0 

1.0 

0.5 

0.0 

17 

-1.0 

0.0 

-0.5 

0.0 

18 

1.0 

0.0 

-0.5 

0.0 

19 

-1.0 

0.0 

0.5 

0.0 

20 

1.0 

0.0 

0.5 

0.0 

21 

0.0 

-1.0 

0.0 

-0.5 

22 

0.0 

1.0 

0.0 

-0.5 

23 

0.0 

-1.0 

0.0 

0.5 

24 

0.0 

1.0 

0.0 

0.5 

25 

0.0 

0.0 

0.0 

0.0 

26 

-0.5 

-1.0 

0.0 

0.0 

27 

0.5 

-1.0 

0.0 

0.0 

28 

-0.5 

1.0 

0.0 

0.0 

29 

0.5 

1.0 

0.0 

0.0 

30 

0.0 

0.0 

-0.5 

-1.0 

31 

0.0 

0.0 

0.5 

-1.0 

32 

0.0 

0.0 

-0.5 

1.0 

33 

0.0 

0.0 

0.5 

1.0 

34 

-0.5 

0.0 

0.0 

-1.0 

35 

0.5 

0.0 

0.0 

-1.0 

36 

-0.5 

0,0 

0.0 

1.0 

37 

0.5 

0.0 

0.0 

1.0 

38 

0.0 

-0.5 

-1.0 

0.0 

39 

0.0 

0.5 

-1.0 

0.0 

40 

0.0 

-0.5 

1.0 

0.0 
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41 

0. 

.0 

0. 

.5 

1 

.0 

0, 

,0 

42 

-0. 

,5 

0. 

,0 

-1 

.0 

0. 

,0 

43 

0. 

,5 

0. 

.0 

-1 

.0 

0, 

,0 

44 

-0. 

.5 

0. 

.0 

1 

.0 

0. 

.0 

45 

0. 

,5 

0. 

.0 

1 

.0 

0. 

.0 

46 

0, 

.0 

-0. 

,5 

0 

.0 

-1, 

,0 

47 

0, 

.0 

0, 

.5 

0 

.0 

-1, 

.0 

48 

0, 

.0 

-0, 

,5 

0 

.0 

1, 

.0 

49 

0, 

.0 

0, 

.5 

0 

.0 

1, 

,0 

50 

0, 

.0 

0, 

.0 

0 

.0 

0, 

,0 

!!!!!!!!!!!!!!!!!!  DECISION  VARIABLES  !!!!!!!!!!!!!!!!!! 
If  OF  DECISION  VARIABLE  SET  CHANGES  :  2 


Decision  var  set  f{  1 
4.0000 
0.0000 
0.6667 
6.6667 

li  of  occurances  of  this  basis:  1003 

%  of  overall  occurance  is:  100.300 

AVE.  OPTIMUM  IS:  215.0843 

Bias  opt  est  (Ave  Opt-True  Opt):  -0.9157 

The  population  variance  is:  3.7599 

The  maximum  value  is:  220.5652 

The  minimum  value  is:  209.3271 

Decision  var  set  H  2 
2.0000 
8.0000 
0.0000 
2.0000 

of  occurances  of  this  basis:  1000 

%  of  overall  occurance  is:  100.000 
AVE.  OPTIMUM  IS:  215.6013 

Bias  opt  est  (Ave  Opt-True  Opt):  -0.3987 

The  population  variance  is:  3.8427 

The  maximum  value  is:  221.0208 

The  minimum  value  is:  209.2253 


!!!!!!!!!!!!!!!!!!!  THE  EXTREME  POINT  !!!!!!!!!!!!!!!!!!! 
NUMBER  OF  EXT  POINTS  VISITED  IS:  2 

Extreme  Point  il  1 
4.0000 
0.0000 
0.6667 
6.6667 

(I  OF  EXT  PT  VISITS  ARE:  1003 

%  OF  OVERALL  IS:  100.300 

AVE.  OPTIMUM  IS:  215.0843 

Bias  opt  est  (Ave  Opt-True  Opt):  -0.9157 
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TRUE  Z*  WITH  TRUE  C  IS:  215.3333 

TRUE  BIAS  (Z*  -  Z  optimal):  -0.6667 

Difference  between  expected  optimal  and  true 
The  population  variance  is:  3.7599 

The  maximum  value  is:  220.5652 

The  minimum  value  is:  209.3271 

Extreme  Point  §  2 

2.0000 
8.0000 
0.0000 
2.0000 

f  OF  EXT  PT  VISITS  ARE:  1000 

%  OF  OVERALL  IS:  100.000 

AVE.  OPTIMUM  IS:  215.6013 

Bias  opt  est  (Ave  Opt-True  Opt):  -0.3987 

TRUE  Z*  WITH  TRUE  C  IS:  216.0000 

TRUE  BIAS  (Z*  -  Z  optimal):  0.0000 

Difference  between  expected  optimal  and  true 
The  population  variance  is:  3.8427 

The  maximum  value  is:  221.0208 

The  minimum  value  is:  209.2253 


!!!!!!!!!!!!!!!!!!!!  OVERALL  RESULTS  !!!!!!!!!! 


Overall  Mean  Optimum  is:  2.8948 

Overall  Bias  (Ave  Opt-True  Opt):  -213.1052 
The  overall  sample  variance  is:  614.7244 

The  overall  maximum  value  is:  221.0208 

The  overall  minimum  vaJ.ue  is:  0.0000 


AVE  OPTIMAL  PER  OBJECTIVE  FUNCTION 
Mean  of  Mean  Opt  per  obj  funct;  2.8948 

The  sample  var(mean  opt):  0.0068 

The  maximum  value  is:  4.3585 

The  minimum  value  is:  2.8164 


0.2490 


0.3987 


M  I  I  I  I  I  I  I  I  I  I 
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NOISE  MULTIPLIER  (SD)  FOR  NORMAL  NOISE:  5.00000 

STANDARD  ERROR  IS  :  1.25000 

This  is  a  screening  run 

Using  a  Design  to  sample 
#  Times  failed  to  sample  true  extreme  pt  0 
Standard  deviation  mult  set  =  1.50000  2.75000 

4.00000 

Number  of  Objective  Function  Samples:  1000 

Number  of  Runs  per  Obj  Function:  149 

Total  Number  of  Points  Tested:  149000 

The  True  Objective  Function: 

const  10.00000+  15.0000  17.0000  18.0000 

20.0000 


Sample 

Generated  Objective  Function 

8  . 

95856 

15.9861 

19.6851 

19.5614 

16 

.8254 

Constraint  Matrix 

1 

2 

3 

4 

1 

1 

1 

2 

1 

2 

2 

1 

-1 

1 

3 

1 

1 

1 

2 

The  RHS  is: 

12.0000 

14.0000 

10.00000 

*True 

Optimal 

Answer : 

216.000 

*True 

optimal 

Extreme  Point: 

2. 

00000 

8.00000 

0.  2.00000 

Design  Matrix 

1 

2 

3 

4 

5 

1 

1 

-1 

-1 

-1 

-1 

2 

1 

1 

-1 

-1 

-1 

3 

1 

-1 

1 

-1 

-1 

4 

1 

1 

1 

-1 

-1 

5 

1 

-1 

-1 

1 

-1 

6 

1 

1 

-1 

1 

-1 

7 

1 

-1 

1 

1 

-1 

8 

1 

1 

1 

1 

-1 

9 

1 

-1 

-1 

-1 

1 

10 

1 

!■ 

-1 

-1 

1 

11 

1 

-1 

1 

-1 

1 

12 

3 

1 

1 

-1 

1 

13 

3 

-1 

-1 

1 

1 

14 

1 

1 

-1 

1 

i 

15 

1 

-1 

1 

1 

1 

16 

1 

1 

1 

1 

1 

Sample 

response  varaible 

Y: 

-70 

.7948 

-33.0653 

-25.8214 

10.42741 

-18  . 

,5379 

17.4076  13 

.7993 

43 

.6501 

-25.0401 

1.01661 

11.4050 

47  . 

,0496 

3.49555  39 

.7061 

55.2738  73.3654 

Sample  Variance-Covariance  Matrix 

12  3  4 
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1 

2.252 

0.000 

0.000 

0.000 

2 

0.000 

2.252 

0.000 

0.000 

3 

0.000 

0.000 

2.252 

0.000 

4 

0.000 

0.000 

0.000 

2.252 

Sample  Cholesky  Factorization  Matrix 

1 

2 

3 

4 

1 

1.501 

0.000 

0.000 

0.000 

2 

1.501 

0.000 

0.000 

3 

1.501 

0.000 

4 

1.501 

SAMPLING 

DESING 

(BOX-BEHNKEN) 

1 

2 

3 

4 

1 

-1.0 

-0.5 

0.0 

0.0 

2 

1.0 

-0.5 

0.0 

0.0 

3 

-1.0 

0.5 

0.0 

0.0 

4 

1.0 

0.5 

0.0 

0.0 

5 

0.0 

0.0 

-1.0 

-0.5 

6 

0.0 

0.0 

1.0 

-0.5 

7 

0.0 

0.0 

-1.0 

0.5 

8 

0.0 

0.0 

1.0 

0.5 

9 

-1.0 

0.0 

0.0 

-0.5 

10 

1.0 

0.0 

0.0 

-0.5 

11 

-1.0 

0.0 

0.0 

0.5 

12 

1.0 

0.0 

0.0 

0.5 

13 

0.0 

-1.0 

-0.5 

0.0 

14 

0.0 

1.0 

-0.5 

0.0 

15 

0.0 

-1.0 

0.5 

0.0 

16 

0.0 

1.0 

0.5 

0.0 

17 

-1.0 

0.0 

-0.5 

0.0 

18 

1.0 

0.0 

-0.5 

0.0 

19 

-1.0 

0.0 

0.5 

0.0 

20 

1.0 

0.0 

0.5 

0.0 

21 

0.0 

-1.0 

0.0 

-0.5 

22 

0.0 

1.0 

0.0 

-0.5 

23 

0.0 

-1.0 

0.0 

0.5 

24 

0.0 

1.0 

0.0 

0.5 

25 

0.0 

0.0 

0.0 

0.0 

26 

-0.5 

-1.0 

0.0 

0.0 

27 

0.5 

-]  .0 

0.0 

0.0 

28 

-0.5 

1.0 

0.0 

0.0 

29 

0.5 

1.0 

0.0 

0.0 

30 

0.0 

0.0 

-0.5 

-1.0 

31 

0.0 

0.0 

0.5 

-1.0 

32 

0.0 

0.0 

-0.5 

1.0 

33 

0.0 

0.0 

0.5 

1.0 

34 

-0.5 

0.0 

0.0 

-1.0 

35 

0.5 

0.0 

0.0 

-1.0 

36 

-0.5 

0.0 

0.0 

1.0 

37 

0.5 

0.0 

0.0 

1.0 

38 

0.0 

-0.5 

-1.0 

0.0 

39 

0.0 

0.5 

-1.0 

0.0 

40 

0.0 

-0.5 

1.0 

0.0 
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41 

0. 

,0 

0. 

5 

1. 

.0 

0. 

,0 

42 

-0. 

,5 

0. 

0 

-1. 

.0 

0. 

,0 

43 

0. 

,5 

0. 

0 

-1. 

.0 

0, 

.0 

44 

-0. 

.5 

0. 

0 

1. 

.0 

0, 

.0 

45 

0, 

,5 

0. 

0 

1. 

,0 

0. 

,0 

46 

0. 

.0 

-0. 

5 

0. 

.0 

-1. 

.0 

47 

0, 

,0 

0. 

5 

0. 

.0 

-1. 

.0 

48 

0. 

.0 

-0. 

5 

0. 

.0 

1. 

.0 

49 

0. 

.0 

0. 

5 

0. 

.0 

1. 

.0 

50 

0. 

.0 

0. 

0 

0. 

,0 

0. 

,0 

!  !  !  ! 

1  1  t  t  1  1  1  1  1  1 

!  !  !  ! 

!  DECISION 

VARIABLES 

1  1  I  t  t  1 

1  1  1  t  1  1  1 

!  ! 

!  OF 

DECISION 

VARIABLE  SET  CHANGES  : 

3 

Decision  var  set  if  1 
4.0000 
0.0000 
0.6667 
6.6667 

ft  of  occurences  of  this  basis:  1002 

%  of  overall  occurence  is:  100.200 

AVE.  OPTIMUM  IS:  212.5927 

Bias  opt  est  (Ave  Opt-True  Opt):  -3.4073 

The  population  variance  is:  91.0682 

The  maximum  value  is:  248.6682 

The  minimum  value  is:  180.1154 

Decision  var  set  ff  2 
2.0000 
10.0000 
0.0000 
0.0000 

#  of  occurences  of  this  basis:  290 

%  of  overall  occurence  is:  29.000 

AVE.  OPTIMUM  IS:  218.1246 

Bias  opt  est  (Ave  Opt-True  Opt):  2.1246 

The  population  variance  is:  95.3511 

The  maximum  value  is:  245.0287 

The  minimum  value  is:  192.3946 

Decision  var  set  ff  3 
2.0000 
8.0000 
0.0000 
2.0000 

ff  of  occurences  of  this  basis:  1000 

%  of  overall  occurence  is:  100.000 

AVE.  OPTIMUM  IS:  217.7167 

Bias  opt  est  (Ave  Opt-True  Opt):  1.7166 

The  population  variance  is:  107.4373 

The  maximum  value  is:  252.3086 

The  minimum  value  is:  185.5855 
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!!!!!!!!!!!!!!!!!!!  THE  EXTREiME  POINT  !!!!!!!! 
NUMBER  OF  EXT  POINTS  VISITED  IS:  4 

Extreme  Point  #  1 

4.0000 
0.0000 
0.6667 
6.6667 

#  OF  EXT  PT  VISITS  ARE:  1002 

%  OF  OVERALL  IS:  100.200 

AVE.  OPTIMUM  IS:  212.5927 

Bias  opt  est  (Ave  Opt-True  Opt):  -3.4073 

TRUE  Z*  WITH  TRUE  C  IS:  215.3333 

TRUE  BIAS  (Z*  -  Z  optimal):  -0.6667 

Difference  between  expected  optimal  and  true 
The  population  variance  is:  91.0682 

The  maximum  value  is:  248.6682 

The  minimum  value  is:  180.1154 

Extreme  Point  #  2 

1.0000 
11.0000 
0.0000 
0.0000 

i\  OF  EXT  PT  VISITS  ARE:  287 

%  OF  OVERALL  IS:  28.700 

AVE.  OPTIMUM  IS:  218.0409 

Bias  opt  est  (Ave  Opt-True  Opt):  2.0409 

TRUE  Z*  WITH  TRUE  C  IS:  210.0000 

TRUE  BIAS  (Z*  -  Z  optimal):  -6.0000 

Difference  between  expected  optimal  and  true 
The  population  variance  is:  95.5884 

The  maximum  value  is:  245.0287 

The  minimum  value  is:  192.3946 

Extreme  Point  tl  3 
2.0000 
8.0000 
0.0000 
2.0000 

(f  OF  EXT  PT  VISITS  ARE:  1000 

%  OF  OVERALL  IS:  100.000 

AVE.  OPTIMUM  IS:  217.7167 

Bias  opt  est  (Ave  Opt-True  Opt):  1.7166 

TRUE  Z*  WITH  TRUE  C  IS:  216.0000 

TRUE  BIAS  (Z*  -  Z  optimal):  0.0000 

Difference  between  expected  optimal  and  true 
The  population  variance  is:  107.4373 

The  maximum  value  is:  252.3086 

The  minimum  value  is:  185.5855 


I  !  !  I  I  I  !  !  I  I  I 


2.7406 


-8.0409 


-1.7167 
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Extreme  Point  #  4 

2.0000 
10.0000 
0.0000 
0.0000 

^  OF  EXT  PT  VISITS  ARE:  3 

%  OF  OVERALL  IS:  0.300 

AVE.  OPTIMUM  IS:  226.1286 

Bias  opt  est  (Ave  Opt-True  Opt):  10.1286 

TRUE  Z*  WITH  TRUE  C  IS:  10.0000 

TRUE  BIAS  (Z*  -  Z  optimal):  -206.0000 

Difference  betv/een  expected  optimal  and  true  -216.1286 

The  population  variance  is:  7.9169 

The  maximum  value  is:  229.8806 

The  minimum  value  is:  223.1049 


!!!!!!!!!!!!!!!!!!!!  OVERALL  RESULTS  !!!!!!!!!!!!!!!!!!!!!! 


Overall  Mean  Optimum  is:  4.4164 

Overall  Bias  (Ave  Opt-True  Opt):  -211.5836 
The  overall  sample  variance  is:  932.7665 

The  overall  maximum  value  is:  252.3086 

The  overall  minimum  value  is:  0.0000 


AVE  OPTIMAL  PER  OBJECTIVE  FUNCTION 
Mean  of  Mean  Opt  per  obj  funct:  3.3154 

The  sample  var(mean  opt):  0.4727 

The  maximum  value  is:  5.9063 

The  minimum  value  is:  2.5426 
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NOISE  MULTIPLIER  (SD)  FOR  NORMAL  NOISE:  9.00000 

STANDARD  ERROR  IS  :  2.25000 

This  is  a  screening  run 

Using  a  Design  to  sample 
~  Times  failed  to  sample  true  extreme  pt  1 
Standard  deviation  mult  set  =  1.10000  2.75000 

4.00000 

%  failures  overall:  0.100 

Number  of  Objective  Function  Samples:  1000 

Number  of  Runs  per  Obj  Function:  149 

Total  Number  of  Points  Tested:  149000 

The  True  Objective  Function: 

const  10.00000+  15.0000  17.0000  18.0000 

20.0000 


Sample  Generated  Objective  Function 

11.3725 

17.5542  15.3801 

15.6525 

22 

.4790 

Constraint  Matrix 

1  2 

3 

4 

1 

1  1 

2 

1 

2 

2  1 

-1 

1 

3 

- 

1  1 

1 

2 

The 

RHS  is: 

12.0000  14.0000 

10.00000 

*True  Optimal 

Answer:  216.000 

*True  optimal 

Extreme  Point: 

2.00000 

8.00000  0.  2.00000 

Design  Matrix 

1 

2  3 

4 

5 

i 

1 

-1  -1 

-2 

-1 

2 

1 

1  -1 

-1 

-1 

3 

1 

-1  1 

-1 

-1 

4 

1 

1  1 

-1 

-1 

5 

1 

-1  -1 

1 

-1 

6 

1 

1  -1 

1 

-1 

7 

1 

-1  1 

1 

-1 

8 

1 

1  1 

1 

-1 

9 

1 

-1  -1 

-1 

1 

10 

1 

1  -1 

-1 

1 

11 

3. 

-1  1 

-1 

1 

12 

1 

1  1 

-1 

1 

13 

1 

-1  -1 

1 

1 

14 

1 

1  -1 

1 

1 

15 

1 

-1  1 

1 

1 

16 

1 

1  1 

1 

1 

Sample  response  varaible  Y: 

-37.9410 

-35.1582  -41.4563 

20.0342 

-29 

.4672 

9.90238  2. 

25696 

22.9770 

-29.1399  17.1932 

17.9181 

54 

.3095 

11. 

6460  60 

.9039 

56.7295 

81.2516 

Sample  Variance-Covariance 

Matrix 
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1 

2 

3 

4 

1 

9.011 

0.000 

0.000 

0.000 

2 

0.000 

9.011 

0.000 

0.000 

3 

0.000 

0.000 

9.011 

0.000 

4 

0.000 

0.000 

0.000 

9.011 

Sample  Cholesky  Factorization  Matrix 

1 

2 

3 

4 

1 

3.002 

0.000 

0.000 

0.000 

2 

3.002 

0.000 

0.000 

3 

3.002 

0.000 

4 

3.002 

SAMPLING 

DESING 

(BOX-BEHNKEN) 

1 

2 

3 

4 

1 

-1.0 

-0.5 

0.0 

0.0 

2 

1.0 

-0.5 

0.0 

0.0 

3 

-1.0 

0.5 

0.0 

0.0 

4 

1.0 

0.5 

0.0 

0.0 

5 

0.0 

0.0 

-1.0 

-0.5 

6 

0.0 

0.0 

1.0 

-0.5 

7 

0.0 

0.0 

-1.0 

0.5 

8 

0.0 

0.0 

1.0 

0.5 

9 

-1.0 

0.0 

0.0 

-0.5 

10 

1.0 

0.0 

0.0 

-0.5 

11 

-1.0 

0.0 

0.0 

0.5 

12 

1.0 

0.0 

0.0 

0.5 

13 

0.0 

-1.0 

-0.5 

0.0 

14 

0.0 

1.0 

-0.5 

0.0 

15 

0.0 

-1.0 

0.5 

0.0 

16 

0.0 

1.0 

0.5 

0.0 

17 

-1.0 

0.0 

-0.5 

0.0 

18 

1.0 

0.0 

-0.5 

0.0 

19 

-1.0 

0.0 

0.5 

0.0 

20 

1.0 

0.0 

0.5 

0.0 

21 

0.0 

-1.0 

0.0 

-0.5 

22 

0.0 

1.0 

0.0 

-0.5 

23 

0.0 

-1.0 

0.0 

0.5 

24 

0.0 

1.0 

0.0 

0.5 

25 

0.0 

0.0 

0.0 

0.0 

26 

-0.5 

-1.0 

0.0 

0.0 

27 

0.5 

-1.0 

0.0 

0.0 

28 

-0.5 

1.0 

0.0 

0.0 

29 

0.5 

1.0 

0.0 

0.0 

30 

0.0 

0.0 

-0.5 

-1.0 

31 

0.0 

0.0 

0.5 

-1.0 

32 

0.0 

0.0 

-0.5 

1.0 

33 

0.0 

0.0 

0.5 

1.0 

34 

-0.5 

0.0 

0.0 

-1.0 

35 

0.5 

0.0 

0.0 

-1.0 

36 

-0.5 

0.0 

0.0 

1.0 

37 

0.5 

0.0 

0.0 

1.0 

38 

0.0 

-0.5 

-1.0 

0.0 

39 

0.0 

0.5 

-1.0 

0.0 
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40 

0. 

,0 

-0.5 

1. 

.0 

0. 

.0 

41 

0. 

,0 

0.5 

1. 

,0 

0. 

,0 

42 

-0. 

.5 

0.0 

-1. 

,0 

0, 

,0 

43 

0, 

,5 

0.0 

-1, 

,0 

0. 

.0 

44 

-0. 

,5 

0.0 

1. 

.0 

0, 

.0 

45 

0. 

,5 

0.0 

1, 

.0 

0; 

,0 

46 

0. 

,0 

-0.5 

0. 

.0 

-1, 

,0 

47 

0. 

.0 

0.5 

0, 

,0 

-1, 

,0 

48 

0, 

.0 

-0.5 

0, 

.0 

1. 

.0 

49 

0, 

.0 

0.5 

0, 

.0 

1. 

.0 

50 

0, 

.0 

0.0 

0, 

.0 

0, 

.0 

1  f  1  1  1  1 

1  M  1  1  1 

MM 

!  DECISION 

VARIABLES 

!  !  !  !  ! 

1  1  M  1  1  1 

!  ! 

#  OP  DECISION  VARIABLE  SET  CHANGES  :  4 


Decision  var  set  1 
4.0000 
0.0000 
0.6667 
6.6667 

#  of  occurances  of  this  basis:  1000 

%  of  overall  occurance  is:  100.000 
AVE.  OPTIMUM  IS:  210.3186 

Bias  opt  est  (Ave  Opt-True  Opt):  -5.6814 

The  population  variance  is:  306.0520 

The  maximum  value  is:  264.7609 

The  minimum  value  is:  158.1760 

Decision  var  set  2 
1.0000 
11.0000 
0.0000 
0.0000 

i(  of  occurances  of  this  basis:  620 

%  of  overall  occurance  is:  62.000 

AVE.  OPTIMUM  IS:  220.3541 

Bias  opt  est  (Ave  Opt-True  Opt):  4.3540 

The  population  variance  is:  334.3401 

The  maximum  value  is:  267.6013 

The  minimum  value  is:  169.1727 

Decision  var  set  #  3 

2.0000 
8.0000 
0.0000 
2.0000 

#  of  occurances  of  this  basis:  999 

%  of  overall  occurance  is:  99.900 
AVE.  OPTIMUM  IS:  220.8882 

Bias  opt  est  (Ave  Opt-True  Opt):  4.8881 

The  population  variance  is:  442.1299 

The  maximum  value  is:  276.3776 

The  minimum  value  is:  161.5952 
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Decision  var  set  #  4 

8.0000 
0.0000 
2.0000 
0.0000 

#  of  occurances  of  this  basis;  19 

%  of  overall  occurance  is:  1.900 

AVE.  OPTIMUM  IS:  208.8170 

Bias  opt  est  (Ave  Opt-True  Opt):  -7.1830 

The  population  variance  is:  380.0007 

The  maximum  value  is:  260.9492 

The  minimum  value  is:  171.0470 


!!!!!!!!!!!!!!!!!!!  THE  EXTREME  POINT  !!!!!!!! 
NUMBER  OF  EXT  POINTS  VISITED  IS:  5 

Extreme  Point  #  1 

4.0000 
0.0000 
0.6667 
6.6667 

jf  OF  EXT  PT  VISITS  ARE:  1000 

%  OP  OVERALL  IS:  100.000 

AVE.  OPTIMUM  IS:  210.3186 

Bias  opt  est  (Ave  Opt-True  Opt):  -5.6814 

TRUE  Z*  WITH  TRUE  C  IS:  215.3333 

TRUE  BIAS  (Z*  -  Z  optimal):  -0.6667 

Difference  between  expected  optimal  and  true 
The  population  variance  is:  306.0520 

The  maximum  value  is:  264.7609 

The  minimum  value  is:  158.1760 

Extreme  Point  #  2 

1.0000 
11.0000 
0.0000 
0.0000 

H  OF  EXT  PT  VISITS  ARE:  558 

%  OF  OVERALL  IS:  55.800 

AVE.  OPTIMUM  IS:  220.3984 

Bias  opt  est  (Ave  Opt-True  Opt):  4.3984 

TRUE  Z*  WITH  TRUE  C  IS:  212.0000 

TRUE  BIAS  (Z*  -  Z  optimal):  -4.0000 

Difference  between  expected  optimal  and  true 
The  population  variance  is:  336.2216 

The  maximum  value  is:  267.6013 

The  minimum  value  is:  169.1727 

Extreme  Point  #  3 

2.0000 


I  I  I  I  I  I  I  I  I  I  I 


5.0147 


-8.3984 
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8.0000 

0.0000 

2.0000 

#  OF  EXT  PT  VISITS  ARE:  999 

%  OF  OVERALL  IS:  99.900- 

AVE.  OPTIMUM  IS:  220.8882 

Bias  opt  est  (Ave  Opt -True  Opt):  4.8881 

TRUE  Z*  WITH  TRUE  C  IS:  216.0000 

TRUE  BIAS  (Z*  -  Z  optimal):  0.0000 

Difference  between  expected  optimal  and  true  -4.8881 

The  population  variance  is:  442.1299 

The  maximum  value  is:  276.3776 

The  minimum  value  is:  161.5952 

Extreme  Point  il  4 
2.0000 
10.0000 
0.0000 
0.0000 

OP  EXT  FT  VISITS  ARE:  62 

%  OF  OVERALL  IS:  6.200 

AVE.  OPTIMUM  IS:  219.9551 

Bias  opt  est  (Ave  Opt-True  Opt):  3.9550 

TRUE  Z*  WITH  TRUE  C  IS:  166.0000 

TRUE  BIAS  (Z*  -  Z  optimal):  -50.0000 

Difference  between  expected  optimal  and  true  -53.9551 

The  population  variance  is:  317.2269 

The  maximum  value  is:  257.6813 

The  minimum  value  is:  174.5772 

Extreme  Point  #  5 

8.0000 
0.0000 
2.0000 
0.0000 

a  OF  EXT  PT  VISITS  ARE:  19 

%  OF  OVERALL  IS:  1.900 

AVE.  OPTIMUM  IS:  208.8170 

Bias  opt  est  (Ave  Opt-True  Opt):  -7.1830 

TRUE  Z*  WITH  TRUE  C  IS:  10.0000 

TRUE  BIAS  (Z*  -  Z  optimal):  -206.0000 

Difference  between  expected  optimal  and  true  -198.8170 

The  population  variance  is:  380.0007 

The  maximum  value  is:  260.9492 

The  minimum  value  is:  171.0470 


!!!!!!!!!!!!!!!!!!!!  OVERALL  RESULTS  !!!!!!!!!!!!!!!!!!!!!! 

Overall  Mean  Optimum  is:  6.1288 

Overall  Bias  (Ave  Opt-True  Opt):  -209.8712 


89 


The  overall  sample  variance  is,: 
The  overall  maximum  value  is : 
The  overall  minimum  value  is: 


AVE  OPTIMAL  PER  OBJECTIVE  FUNCTION 
Mean  of  Mean  Opt  per  obj  funct: 

The  sample  var(mean  opt):  0. 
The  maxim'um  value  is:  8.1238 
The  minimum  value  is:  2.3285 


1298.7844 

276.3776 

0.0000 


3.8361 

8616 
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NOISE  MULTIPLIER  (SD)  FOR  NORMAL  NOISE:  13.0000 

STANDARD  ERROR  IS  :  3.25000 

This  is  a  screening  run 

Using  a  Design  to  sample 
#  Times  failed  to  sample  true  extreme  pt  5 
Standard  deviation  mult  set  =  1.50000  2.75000 

4.00000 

%  failures  overall:  0.500 

Number  of  Objective  Function  Samples:  1000 

Number  of  Runs  per  Obj  Function:  149 


Total 

Number  of  Points  Tested:  149000 

The  True  Objective  Function: 

const 

10.00000+  15. 

0000  17. 

0000  18. 

,0000 

20.0000 

Sample  Generated  Objective  Function 

11 

.8447  15.9522 

17.5825 

18.9583 

16 

.7721 

Constraint  Matrix 

1 

2 

3 

4 

1 

1 

1 

2 

1 

2 

2 

1 

-1 

1 

3 

-1 

1 

1 

2 

The  RHS  is:  12.0000 

14.0000 

10.00000 

*True 

Optimal  Answer: 

216.000 

*True 

optimal  Extreme  Point: 

2. 

00000  8.00000 

0.  2.00000 

Design  Matrix 

1  2 

3 

4 

5 

1 

1  -1 

-1 

-1 

-1 

2 

1  1 

-1 

-1 

-1 

3 

1  -1 

1 

-1 

-1 

4 

1  1 

1 

-1 

-1 

5 

1  -1 

-1 

1 

-1 

6 

1  1 

-1 

1 

-1 

7 

1  -1 

1 

1 

-1 

8 

1  1 

1 

1 

-1 

9 

1  -1 

-1 

-1 

1 

10 

1  1 

-1 

-1 

1 

11 

1  -1 

1 

-1 

1 

12 

1  1 

1 

-1 

1 

13 

1  -1 

-1 

1 

1 

14 

1  1 

-1 

1 

1 

15 

1  -1 

1 

1 

1 

16 

1  1 

1 

1 

1 

Sample 

response  varaible 

Y: 

-45 

.8560  -25.5360 

-27.8729 

16.8729 

-22, 

,9366 

2.03671  12.6677 

51 

.2050  -26.2992 

4.34243 

-8.73610 

56, 

.1767 

25.7303  42.6163 

60.4428  74.6617 

Sample  Variance-Covariance  Matrix 
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1 

2 

3 

4 

1 

7.179 

0.000 

0.000 

0.000 

2 

0.000 

7.179 

0.000  • 

0.000 

3 

0.000 

0.000 

7.179 

0.000 

4 

0.000 

0.000 

0.000 

7.179 

Sample  Cholesky  Factorization  Matrix 

1 

2 

3 

4 

1 

2.679 

0.000 

0.000 

0.000 

2 

2.679 

0.000 

0.000 

3 

2.679 

0.000 

4 

2.679 

SAMPLING 

DESING 

(BOX-BEHNKEN) 

1 

2 

3 

4 

1 

-1.0 

-0.5 

0.0 

0.0 

2 

1.0 

-0.5 

0.0 

0.0 

3 

-1.0 

0.5 

0.0 

0.0 

4 

1.0 

0.5 

0.0 

0.0 

5 

0.0 

0.0 

-1.0 

-0.5 

6 

0.0 

0.0 

1.0 

-0.5 

7 

0.0 

0.0 

-1.0 

0.5 

8 

0.0 

0.0 

1.0 

0.5 

9 

-1.0 

0.0 

0.0 

-0.5 

10 

1.0 

0.0 

0.0 

-0.5 

11 

-1.0 

0.0 

0.0 

0.5 

12 

1.0 

0.0 

0.0 

0.5 

13 

0.0 

-1.0 

-0.5 

0.0 

14 

0.0 

1.0 

-0.5 

0.0 

15 

0.0 

-1.0 

0.5 

0.0 

16 

0.0 

1.0 

0.5 

0.0 

17 

-1.0 

0.0 

-0.5 

0.0 

18 

1.0 

0.0 

-0.5 

0.0 

19 

-1.0 

0.0 

0.5 

0.0 

20 

1.0 

0.0 

0.5 

0.0 

21 

0.0 

-1.0 

0.0 

-0.5 

22 

0.0 

1.0 

0.0 

-0.5 

23 

0.0 

-1.0 

0.0 

0.5 

24 

0.0 

1.0 

0.0 

0.5 

25 

0.0 

0.0 

0.0 

0.0 

26 

-0.5 

-1.0 

0.0 

0.0 

27 

0.5 

-1.0 

0.0 

0.0 

28 

-0.5 

1.0 

0.0 

0.0 

29 

0.5 

1.0 

0.0 

0.0 

30 

0.0 

0.0 

-0.5 

-1.0 

31 

0.0 

0.0 

0.5 

-1.0 

32 

0.0 

0.0 

-0.5 

1.0 

33 

0.0 

0.0 

0.5 

1.0 

34 

-0.5 

0.0 

0.0 

-1.0 

35 

0.5 

0.0 

0.0 

-1.0 

36 

-0.5 

0.0 

0.0 

1.0 

37 

0.5 

0.0 

0.0 

1.0 

38 

0.0 

-0.5 

-1.0 

0.0 

39 

0.0 

0.5 

-1.0 

0.0 
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4  a- 

0.0 

-0.5 

■i.a 

0.0 

41 

0.0 

0.5 

1.0 

0.0 

42 

-6.5 

0.0 

’  -1.0 

0..0 

43 

0.5 

0.0 

-1.0 

0.0 

44 

-0.5 

0.0 

1.0 

0.0 

45 

0.5 

0.0 

1.0 

0.0 

46 

O.d 

-0.5 

0.0 

-1.0 

47 

0.0 

0;5 

0.0 

-1.0 

48 

0.0 

-0.5 

0.0 

1.0 

49 

0.0 

0.5 

0.0' 

1.0 

50 

0.0 

0.0 

6.0 

0.0 

1  1  1  1  1  1  1 

1  1  I  1  1  1. 1  1  1 

!  DECISION 

VARIABLES  ! ! ! ! ! 

1  1  1  1  1  I  1  1  i 

#  OF  DECISION  VARIABLE  SET  CHANGES  ;  6 


Decision  var  set  #  1 

4.0000 
0.0000 
0.6667 
6.6667 

#  of  occurances  of  this  basis:  999 

%  of  overall  occurance  is:  99.900 
AVE.  OPTIMUM  IS:  208.5842 

Bias  opt  est  (Ave  Opt-True  Opt):  -7.4158 

The  population  variance  is:  651.1973 

The  maximum  value  is:  288.6066 

The  minimum  value  is:  137.1993 

Decision  var  set  if  2 
2.0000 
8.0000 
0.0000 
2.0000 

ft  of  occurances  of  this  basis:  995 

%  of  overall  occurance  is:  99.500 

AVE.  OPTIMUM  IS:  222.7262 

Bias  opt  est  (Ave  Opt-True  Opt):  6.7262 

The  population  variance  is:  841.8345 

The  maximum  value  is:  310.9501 

The  minimum  value  is:  145.5584 

Decision  var  set  #  3 

1.0000 
11.0000 
0.0000 
0.0000 

#  of  occurances  of  this  basis:  806 

%  of  overall  occurance  is:  80.600 
AVE.  OPTIMUM  IS:  223.1638 

Bias  opt  est  (Ave  Opt-True  Opt):  7.1638 

The  population  variance  is:  721.3096 

The  maximum  value  is:  308.3556 

The  minimum  value  is:  150.9656 
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Decision  var  set  #  4 

8.0000 
0.0000 
2.0000 
0.0000 

#  of  occurances  of  this  basis:  103 

%  of  overall  occurance  is:  10.300 

AVE.  OPTIMUM  IS:  225  2866 

Bias  opt  est  (Ave  Opt-True  Opt):  9.2865 

The  population  variance  is:  760.8097 

The  maximum  value  is:  283.2927 

The  minimum  value  is:  151.0016 

Decision  var  set  #  5 

0.0000 
0.0000 

4.6667 

2.6667 

of  occurances  of  this  basis:  7 

%  of  overall  occurance  is:  0.700 

AVE.  OPTIMUM  IS:  168.4372 

Bias  opt  est  (Ave  Opt-True  Opt):  -47.5629 

The  population  variance  is:  661.1362 

The  maximum  value  is:  219.5729 

The  minimum  value  is:  139.5171 

Decision  var  set  #  6 

0.0000 
8.0000 
2.0000 
0.0000 

^  of  occurances  of  this  basis:  5 

%  of  overall  occurance  is:  0.500 

AVE.  OPTIMUM  IS:  163.4464 

Bias  opt  est  (Ave  Opt-True  Opt):  -52.5536 

The  population  variance  is:  165.5492 

The  maximum  value  is:  187.7662 

The  minimum  value  is:  149.7473 


!!!!!!!!!!!!!!!!!!!  THE  EXTREME  POINT  !!!!!!!!!!!!!!!!!!! 
NUMBER  OF  EXT  POINTS  VISITED  IS:  7 

Extreme  Point  #  1 

4.0000 
0.0000 
0.6667 

6.6667 

#  OF  EXT  PT  VISITS  ARE:  999 

%  OF  OVERALL  IS:  99.900 

AVE.  OPTIMUM  IS:  208.5842 

Bias  opt  est  (Ave  Opt-True  Opt):  -7.4158 
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TRUE  Z*  WITH  TRUE  C  IS*.  215.3333 

TRUE  BIAS  (Z*  -  Z  optimal):  -0.6667 

Difference  between  expected  optimal  and  true  6.7491 

The  population  variance  is:  651.1973 

The  maximum  value  is:  288.6066 

The  minimum  value  is:  137.1993 

Extreme  Point  #  2 

2.0000 
8. .0000 
0.0000 
2.0000 

#  OF  EXT  PT  VISITS  ARE:  995 

%  OF  OVERALL  IS:  99.500 

AVE.  OPTIMUM  IS:  222.7262 

Bias  opt  est  (Ave  Opt-True  Opt):  6.7262 

TRUE  Z*  WITH  TRUE  C  IS:  216.0000 

TRUE  BIAS  (Z*  -  Z  optimal):  0.0000 

Difference  between  expected  optimal  and  true  -6.7262 

The  population  variance  is:  841.8345 

The  maximum  value  is:  310.9501 

The  minimum  value  is:  145.5584 

Extreme  Point  #  3 

1.0000 
11.0000 
0.0000 
0.0000 

a  OF  EXT  PT  VISITS  ARE:  661 

%  OP  OVERALL  IS:  66.100 

AVE.  OPTIMUM  IS:  221.9296 

Bias  opt  est  (Ave  Opt-True  Opt):  5.9295 

TRUE  Z*  WITH  TRUE  C  IS:  212.0000 

TRUE  BIAS  (Z*  -  Z  optimal):  -4.0000 

Difference  between  expected  optimal  and  true  -9.9296 

The  population  variance  is:  700.3993 

The  maximum  value  is:  297.5756 

The  minimum  value  is:  153.1106 

Extreme  Point  #  4 

2.0000 
10.0000 
0.0000 
0.0000 

if  OF  EXT  PT  VISITS  ARE:  145 

%  OF  OVERALL  IS:  14.500 

AVE.  OPTIMUM  IS:  228.7902 

Bias  opt  est  (Ave  Opt-True  Opt):  12.7901 

TRUE  Z*  WITH  TRUE  C  IS:  166.0000 

TRUE  BIAS  (Z*  -  Z  optimal):  -50.0000 

Difference  between  expected  optimal  and  true  -62.7902 
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The  population  variance  is:  778.0309 

The  maximum  value  is:  308.3556 

The  minimum  value  is:  150.9656 

Extreme  Point  #  5 

8.0000 
0.0000 
2.0000 
0.0000 

#  OF  EXT  PT  VISITS  ARE:  103 

%  OF  OVERALL  IS:  10.300 
AVE.  OPTIMUM  IS:  225.2866 

Bias  opt  est  (Ave  Opt-True  Opt):  9.2865 

TRUE  Z*  WITH  TRUE  C  IS:  147.3333 

TRUE  BIAS  (Z*  -  Z  optimal):  -68.6667 

Difference  between  expected  optimal  and  true  -77.9532 

The  population  variance  is:  760.8097 

The  maximum  value  is:  283.2927 

The  minimum  value  is:  151.0016 

Extreme  Point  j|  6 
0.0000 
0.0000 

4.6667 

2.6667 

#  OF  EXT  PT  VISITS  ARE:  7 

%  OF  OVERALL  IS:  0.700 

AVE.  OPTIMUM  IS:  168.4372 

Bias  opt  est  (Ave  Opt-True  Opt):  -47.5629 

TRUE  Z*  WITH  TRUE  C  IS:  182.0000 

TRUE  BIAS  (Z*  -  Z  optimal):  -34.0000 

Difference  between  expected  optimal  and  true  13.5628 

The  population  variance  is:  661.1362 

The  maximum  value  is:  219.5729 

The  minimum  value  is:  139.5171 

Extreme  Point  (f  7 
0.0000 
8.0000 
2.0000 
0.0000 

#  OF  EXT  PT  VISITS  ARE:  5 

%  OF  OVERALL  IS:  0.500 

AVE.  OPTIMUM  IS:  163.4464 

Bias  opt  est  (Ave  Opt-True  Opt):  -52.5536 

TRUE  Z*  WITH  TRUE  C  IS:  10.0000 

TRUE  BIAS  (Z*  -  Z  optimal):  -206.0000 

Difference  between  expected  optimal  and  true  -153.4464 

The  population  variance  is:  165.5492 

The  maximum  value  is:  187.7662 

The  minimum  value  is:  149.7473 
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!!!!!!!!!!!!!!!!!!!!  OVERALL  RESULTS  !!!!<!!!!!!!!!!!!!!!!! 


Overall  Mean  Optimum  is:  7.9645 

Overall  Bias  (Ave  Opt-True  Opt):  -208.0355 
The  overall  sample  variance  is:  1695.8068 

The  overall  maximum  value  is:  310.9501 

The  overall  minimum  value  is:  0.0000 


AVE  OPTIMAL  PER  OBJECTIVE  FUNCTION 
Mean  of  Mean  Opt  per  obj  funct:  4.2621 

The  sample  var(mean  opt):  1.5189 

The  maximum  value  is:  9.1224 

The  minimum  value  is:  2.0629 
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NOISE  MULTIPLIER  (SD)  FOR  NORMAL  NOISE:  17.0000 

STANDARD  ERROR  IS  :  4.25000 

This  is  a  screening  run 

Using  a  Design  to  sample 
#  Times  failed  to  sample  true  extreme  pt  19 
Standard  deviation  mult  set  =  1.50000  2.75000 

4.00000 

%  failures  overall:  1.900 

Number  of  Objective  Function  Samples:  1000 

Number  of  Runs  per  Obj  Function:  149 

Total  Number  of  Points  Tested:  149000 

The  True  Objective  Function: 

const  10.00000+  15.0000  17.0000  18.0000 

20.0000 

Sample  Generated  Objective  Function 


8.30373 

10.6681 

21.0220 

19.5529 

19 

Constraint  Matrix 

1 

2 

3 

4 

1 

1 

1 

2 

1 

2 

2 

1 

-1 

1 

3 

1 

1 

1 

2 

The  RHS  is: 

12.0000 

14.0000 

10.00000 

*True  Optimal  i 

^nswer : 

216.000 

*True  optimal  Extreme  Point: 


2.00000  8.00000  0.  2.00000 


Design  Matrix 

1 

2 

3 

4 

5 

1 

1 

-1 

-1 

-1 

-1 

2 

1 

1 

-1 

-1 

-1 

3 

1 

-1 

1 

-1 

-1 

4 

1 

1 

1 

-1 

-1 

5 

1 

-1 

-1 

1 

-1 

6 

1 

1 

-1 

1 

-1 

7 

1 

-1 

1 

1 

-1 

8 

1 

1 

1 

1 

-1 

9 

1 

-1 

-1 

-1 

1 

10 

1 

1 

-1 

-1 

1 

11 

1 

-1 

1 

-1 

1 

12 

1 

1 

1 

-1 

1 

13 

1 

-1 

-1 

1 

1 

14 

1 

1 

-1 

1 

1 

15 

1 

-1 

1 

1 

1 

16 

1 

1 

1 

1 

1 

Sample  response  varaible 

Y: 

-54.3485 

-48.7402 

-26.4912 

16.6902 

-35 

.1069 

4.32857  2. 

29300 

51.6007 

-26.1163 

-5.40654 

16.1785 

38 

.2412 

32. 

9653  30 

.6787 

71.7110 

64.3823 

Sample  Variance-Covariance  Matrix 
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1 

2 

3 

4 

1 

10.04 

0.00 

0.00 

0.00 

2 

0.00 

10.04 

0.00 

0.00 

3 

0.00 

0.00 

10.04 

0.00 

4 

0.00 

0.00 

0.00 

10.04 

Sample  Cholesky  Factorization  Matrix 
12  3 

4 

1 

3.169 

0.000 

0.000 

0.000 

2 

3.169 

0.000 

0.000 

3 

4 

SAMPLING 

1 

DESING 

2 

3.169 

(BOX-BEHNKEN) 

3 

0.000 

3.169 

4 

1 

-1.0 

-0.5 

0.0 

0.0 

2 

1.0 

-0.5 

0.0 

0.0 

3 

-1.0 

0.5 

0.0 

0.0 

4 

1.0 

0.5 

0.0 

0.0 

5 

0.0 

0.0 

-1.0 

-0.5 

6 

0.0 

0.0 

1.0 

-0.5 

7 

0.0 

0.0 

-1.0 

0.5 

8 

0.0 

0.0 

1.0 

0.5 

9 

-1.0 

0.0 

0.0 

-0.5 

10 

1.0 

0.0 

0.0 

-0.5 

11 

-1.0 

0.0 

0.0 

0.5 

12 

1.0 

0.0 

0.0 

0.5 

13 

0.0 

-1.0 

-0.5 

0.0 

14 

0.0 

1.0 

-0.5 

0.0 

15 

0.0 

-1.0 

0.5 

0.0 

16 

0.0 

1.0 

0.5 

0.0 

17 

-1.0 

0.0 

-0.5 

0.0 

18 

1.0 

0.0 

-0.5 

0.0 

19 

-1.0 

0.0 

0.5 

0.0 

20 

1.0 

0.0 

0.5 

0.0 

21 

0.0 

-1.0 

0.0 

-0.5 

22 

0.0 

1.0 

0.0 

-0.5 

23 

0.0 

-1.0 

0.0 

0.5 

24 

0.0 

1.0 

0.0 

0.5 

25 

0.0 

0.0 

0.0 

0.0 

26 

-0.5 

-1.0 

0.0 

0.0 

27 

0.5 

-1.0 

0.0 

0.0 

28 

-0.5 

1.0 

0.0 

0.0 

29 

0.5 

1.0 

0.0 

0.0 

30 

0.0 

0.0 

-0.5 

-1.0 

31 

0.0 

0.0 

0.5 

-1.0 

32 

0.0 

0.0 

-0.5 

1.0 

33 

0.0 

0.0 

0.5 

1.0 

34 

-0.5 

0.0 

0.0 

-1.0 

35 

0.5 

0.0 

0.0 

-1.0 

36 

-0.5 

0.0 

0.0 

1.0 

37 

0.5 

0.0 

0.0 

1.0 

38 

0.0 

-0.5 

-1.0 

0.0 

39 

0.0 

0.5 

-1.0 

0.0 
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40 

0.0 

-0.5 

1.0 

0.0 

41 

0.0 

0.5 

1.0 

0.0 

42 

-0.5 

0.0 

-1.0 

0.0 

43 

0.5 

0.0 

-1.0 

0.0 

44 

-0.5 

0.0 

1.0 

0.0 

45 

0.5 

0.0 

1.0 

0.0 

46 

0.0 

-0.5 

0.0 

-1.0 

47 

0.0 

0.5 

0.0 

-1.0 

48 

0.0 

-0.5 

0.0 

1.0 

49 

0.0 

0.5 

0.0 

1.0 

50 

0.0 

0.0 

0.0 

0.0 

!!!!!!!!!!!!!!!!!!  DECISION  VARIABLES  !!!!!!!!!!!!!!!!!! 
OF  DECISION  VARIABLE  SET  CHANGES  :  7 


Decision  var  set  I!  1 
4.0000 
0.0000 
0.6667 
6.6667 

of  occurances  of  this  basis:  1001 

%  of  overall  occurance  is:  100.100 
AVE.  OPTIMUM  IS:  207.5924 

Bias  opt  est  (Ave  Opt-True  Opt):  -8.4077 

The  population  variance  is:  973.1924 

The  maximum  value  is:  304.9121 

The  minimum  value  is:  103.5766 

Decision  var  set  #  2 

1.0000 
11.0000 
0.0000 
0.0000 

If  of  occurances  of  this  basis:  912 

%  of  overall  occurance  is:  91.200 

AVE.  OPTIMUM  IS:  225.7134 

Bias  opt  est  (Ave  Opt-True  Opt):  9.7133 

The  population  variance  is:  1225.7222 

The  maximum  value  is:  336.5096 

The  minimum  value  is:  129.2516 

Decision  var  set  #  3 

2.0000 
8.0000 
0.0000 
2.0000 

If  of  occurances  of  this  basis:  981 

%  of  overall  occurance  is:  98.100 

AVE.  OPTIMUM  IS:  227.2527 

Bias  opt  est  (Ave  Opt-True  Opt):  11.2526 

The  population  variance  is:  1432.8307 

The  maximum  value  is:  371.2051 

The  minimum  value  is:  120.2687 
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Decision  var  set  if  4 
8.0000 
0.0000 
2.0000 
0.0000 

#  of  occurances  of  this  basis:  188 

%  of  overall  occurance  is:  18.800 

AVE.  OPTIMUM  IS:  232.0063 

Bias  opt  est  (Ave  Opt-True  Opt):  16.0063 

The  population  variance  is;  1289.7496 

The  maximum  value  is:  319.6638 

The  minimum  value  is:  122.2305 

Decision  var  set  if  5 
0.0000 
8.0000 
2.0000 
0.0000 

If  of  occurances  of  this  basis:  51 

%  of  overall  occurance  is:  5.100 

AVE.  OPTIMUM  IS:  174.5180 

Bias  opt  est  (Ave  Opt-True  Opt):  -41.4821 

The  population  variance  is:  595.1209 

The  maximum  value  is:  222.3652 

The  minimum  value  is:  113.8435 

Decision  var  set  If  6 
0.0000 
0.0000 

4.6667 

2.6667 

If  of  occurances  of  this  basis:  47 

%  of  overall  occurance  is:  4.700 

AVE.  OPTIMUM  IS:  164.8198 

Bias  opt  est  (Ave  Opt-True  Opt):  -51.1802 

The  population  variance  is:  478.3920 

The  maximum  value  is:  223.9326 

The  minimum  value  is:  128.0441 

Decision  var  set  If  7 
0.0000 
0.0000 
6.0000 
0.0000 

a  of  occurances  of  this  basis:  19 

%  of  overall  occurance  is:  1.900 

AVE.  OPTIMUM  IS:  175.7839 

Bias  opt  est  (Ave  Opt-True  Opt):  -40.2161 

The  population  variance  is:  1132.7343 

The  maximum  value  is:  243.3237 

The  minimum  value  is:  127.3116 
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!!!!!!!!!!!!!!!!!!!  THE  EXTREME  POINT  !!!!!!!!!!!!!!!!!!! 
NUMBER  OF  EXT  POINTS  VISITED  IS:  8 

Extreme  Point  #  1 

4.0000 
0.0000 
0.6667 
6.6667 

#  OP  EXT  PT  VISITS  ARE:  1001 

%  OF  OVERALL  IS:  100.100 

AVE.  OPTIMUM  IS:  207.5924 

Bias  opt  est  (Ave  Opt-True  Opt):  -8.4077 

TRUE  Z*  WITH  TRUE  C  IS:  215.3333 

TRUE  BIAS  (Z*  -  Z  optimal):  -0.6667 

Difference  between  expected  optimal  and  true  7.7410 

The  population  variance  is:  973.1924 

The  maximum  value  is:  304.9121 

The  minimum  value  is:  103.5766 

Extreme  Point  if  2 
1.0000 
11.0000 
0.0000 
0.0000 

if  OP  EXT  PT  VISITS  ARE:  736 

%  OF  OVERALL  IS:  73.600 

AVE.  OPTIMUM  IS:  224.2276 

Bias  opt  est  (Ave  Opt-True  Opt):  8.2276 

TRUE  Z*  WITH  TRUE  C  IS:  212.0000 

TRUE  BIAS  (Z*  -  Z  optimal):  -4.0000 

Difference  between  expected  optimal  and  true  -12.2276 

The  population  variance  is:  1152.1857 

The  maximum  value  is:  314.2408 

The  minimum  value  is:  129.2516 

Extreme  Point  if  3 
2.0000 
10.0000 
0.0000 
0.0000 

if  OF  EXT  PT  VISITS  ARE:  176 

%  OF  OVERALL  IS:  17.600 

AVE.  OPTIMUM  IS:  231.9267 

Bias  opt  est  (Ave  Opt-True  Opt):  15.9266 

TRUE  Z*  WITH  TRUE  C  IS:  216.0000 

TRUE  BIAS  (Z*  -  Z  optimal):  0.0000 

Difference  between  expected  optimal  and  true  -15.9266 

The  population  variance  is:  1485.4053 

The  maximum  value  is:  336.5096 

The  minimum  value  is:  129.5463 
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Extreme  Point  4 
2.0000 
8.0000 
0.0000 
2.0000 

OF  EXT  PT  VISITS  ARE:  981 

%  OF  OVERALL  IS:  98.100 
AVE.  OPTIMUM  IS:  227.2527 

Bias  opt  est  (Ave  Opt-True  Opt):  11.2526 

TRUE  Z*  WITH  TRUE  C  IS:  166.0000 

TRUE  BIAS  (Z*  -  Z  optimal):  -50.0000 

Difference  between  expected  optimal  and  true  -61.2527 

The  population  variance  is:  1432.8307 

The  maximum  value  is:  371.2051 

The  minimum  value  is:  120.2687 

Extreme  Point  ||  5 

8.0000 
0.0000 
2.0000 
0.0000 

|{  OP  EXT  PT  VISITS  ARE:  188 

%  OF  OVERALL  IS:  18.800 

AVE.  OPTIMUM  IS:  232.0063 

Bias  opt  est  (Ave  Opt-True  Opt):  16.0063 

TRUE  Z*  WITH  TRUE  C  IS:  182.0000 

TRUE  BIAS  (Z*  -  Z  optimal):  -34.0000 

Difference  between  expected  optimal  and  true  -50.0063 

The  population  variance  is:  1289.7496 

The  maximum  value  is:  319.6638 

The  minimum  value  is:  122.2305 

Extreme  Point  #  6 

0.0000 
8.0000 
2.0000 
0.0000 

H  OF  EXT  PT  VISITS  ARE:  51 

%  OF  OVERALL  IS:  5.100 

AVE.  OPTIMUM  IS:  174.5180 

Bias  opt  est  (Ave  Opt-True  Opt):  -41.4821 

TRUE  Z*  WITH  TRUE  C  IS:  147.3333 

TRUE  BIAS  (Z*  -  Z  optimal):  -68.6667 

Difference  between  expected  optimal  and  true  -27.1846 

The  population  variance  is:  595.1209 

The  maximum  value  is:  222.3652 

The  minimum  value  is:  113.8435 

Extreme  Point  #  7 

0.0000 
0.0000 
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4.6667 

2.6667 

#  OF  EXT  PT  VISITS  ARE:  47 

%  OF  OVERALL  IS:  4.700 

AVE.  OPTIMUM  IS:  164.8198 

Bias  opt  est  (Ave  Opt-True  Opt):  -51.1802 

TRUE  Z*  WITH  TRUE  C  IS:  118.0000 

TRUE  BIAS  (Z*  -  Z  optimal):  -98.0000 

Difference  between  expected  optimal  and  true  -46.8198 

The  population  variance  is:  478.3920 

The  maximum  value  is:  223.9326 

The  minimum  value  is:  128.0441 

Extreme  Point  if  8 
0.0000 
0.0000 
6.0000 
0.0000 

a  OF  EXT  PT  VISITS  ARE:  19 

%  OP  OVERALL  IS:  1.900 

AVE.  OPTIMUM  IS:  175.7839 

Bias  opt  est  (Ave  Opt-True  Opt):  -40.2161 

TRUE  Z*  WITH  TRUE  C  IS:  10.0000 

TRUE  BIAS  (Z*  -  Z  optimal):  -206.0000 

Difference  between  expected  optimal  and  true  -165.7839 

The  population  variance  is:  1132.7343 

The  maximum  value  is:  243.3237 

The  minimum  value  is:  127.3116 


!!!!!!!!!!!!!!!!!!!!  OVERALL  RESULTS  !!!!!!!!!!!!!!!!!!!!!! 


Overall  Mean  Optimum  is:  9.9972 

Overall  Bias  (Ave  Opt-True  Opt):  -206.0028 
The  overall  sample  variance  is:  2134.0967 

The  overall  maximum  value  is:  371.2051 

The  overall  minimum  value  is:  0.0000 


AVE  OPTIMAL  PER  OBJECTIVE  FUNCTION 
Mean  of  Mean  Opt  per  obj  funct:  4.6993 

The  sample  var(mean  opt):  2.0915 

The  maximum  value  is:  10.7572 

The  minimum  value  is:  1.6353 
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